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FOREWORD 


This report, prepared by the Martin Marietta Corporation, 
Denver Division, under Contract NAS8-30690, presents the re- 
sults of an analytical and experimental study of transient 
liquid motion similar to that encountered in orbiting space- 
craft. The study was performed from March 1974 to February 
1975 and was administered by the National Aeronautics and 
Space Administration, George C. Marshall Space Flight Center, 
Huntsville, Alabama, under the direction of Mr. Frank Bugg. 
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ABSTRACT 


This report presents the results of a twofold study of 
transient liquid motion such as that which will be experienced 
during orbital maneuvers by Space Tug. A test program was con- 
ducted in the Martin Marietta Low-g Test Facility involving 
twenty-two drops. Biaxial, low-g accelerations were applied 
to an instrumented, model propellant tank during free-fall 
testing, and forces exerted during liquid reorientation were 
measured and recorded. Photographic records of the liquid re- 
orientation were also made. The test data was used to verify 
a mechanical analog which portrays the liquid as a point mass 
moving on an ellipsoidal constraint surface. The mechanical 
analog was coded into a Fortran IV digital computer program: 
LAMPS, Large AMP litude _Slosh. Test/analytical correlation in- 
dicates that the mechanical analog is capable of predicting the 
overall force trends measured during testing. More work is 
needed, however, to fine-tune the model through better under- 
standing of viscous dissipative forces and improvements to 
liquid motion constraints. 
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I. INTRODUCTION 


In the immediate future, vehicles similar to Space Tug 
will perform orbital maneuvers while carrying a large mass of 
propellant. An in-depth understanding of the interaction forces 
between the propellant and space vehicle is required to properly 
assess the dynamics of these maneuvers, in particular, the dock- 
ing maneuver. During orbital maneuvers, the propellant mass is 
subjected to small accelerations which can induce large ampli- 
tude slosh motion. Due to the relatively large mass of propel- 
lant, the forces exerted on the spacecraft by the moving propel- 
lant may have a significant effect on gross vehicle motion. 
Knowledge of these interaction forces is imperative in the de- 
sign of vehicle orbital control systems and docking mechanisms. 

A two-fold study has been conducted to develop a mechani- 
cal analog to simulate large amplitude liquid motion in a 
container for low-g environments. The primary purpose of the 
model is to simulate the interaction forces between liquids 
and the spacecraft due to moving propellant. The study con- 
sisted of both experimental and analytical tasks. 

The Martin Marietta (Denver) Drop Test Facility was used 
in the conduct of the experimental phase. A test module capa- 
ble of measuring large amplitude slosh forces on a scale model, 
axisymmetric tank was constructed. The test and liquid were 
scaled using dimensional analysis techniques to ensure simula- 
tion of motion representative of a full scale liquid oxygen 
tank. During testing, the module was dropped in the free-fall 
tower (simulating zero G) and small biaxial accelerations were 
applied. The ensuing liquid motion was photographed and two 
dimensional forces were measured and recorded. A total of 
twenty- two tests were conducted. Various tank fill volumes, 
tank orientations and acceleration magnitudes were investigated. 
The test time (Ri2 seconds) corresponds to approximately 15 
seconds of liquid motion in a full size tank. Chapter II de- 
tails the test program. 

A mechanical analog was developed in the analytical phase 
to simulate the observed large amplitude slosh. The analog 
portrays the liquid as a point mass moving on a constraint 
surface which is represented by piecewise continuous elliptical 
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segments. The constraint surface is the locus of liquid center 
of mass locations prescribed by slowly rotating the tank in a 
one-g field. The mechanical analog was implemented in a com- 
puter program, LAMPS (Large AMP litude _Slosh) , which predicts 
force time histories on the tank due to the liquid motion. 
Chapter III presents a derivation of the equations of liquid 
motion and discusses the computer solution. A users guide to 
the program, LAMPS, is also presented. 

In Chapter IV, a comparison between the acquired test data 
and computer predictions for the test configurations is pre- 
sented . 



II- 1 


II . EXPERIMENTAL INVESTIGATION 


The primary objective of the experimental investigation 
was to generate data for correlation with the computer model 
developed under the analytical task. It was, therefore, re- 
quired that the tests simulate the type of propellant motion 
that could occur during docking maneuvers of a spacecraft. In 
addition, it was necessary to measure the forces produced by 
the motion of the liquid. Scaling was used to relate the test 
conditions to a full-size tank. The drop tower, low-g test 
facility was selected as the means of performing the tests and 
a series of 22 drop tests was accomplished. Since the test con- 
ditions used have not been experimentally simulated before, the 
tests add to the basic understanding of the reorientation of 
propellant within a tank. 


A. TEST MATRIX 

The tests simulated the full-size liquid oxygen tank shown 
in Figure II- 1. It has hemispherical end domes and a short cyl- 
indrical barrel section. The following test parameters were 
varied so that each test simulated a somewhat different condi- 
tion under which the reorientation of the liquid occurs: 

1. liquid volume; 

2. inclination of the tank with respect to a vertical 
reference axis, given by the angle 9X in Figure II- 1; 

3. magnitude of the acceleration acting on the tank: 
the acceleration is composed of two components, an 
axial acceleration and a lateral acceleration. 

A lateral acceleration was applied in order to force the 
liquid to follow the tank wall. This avoids producing insta- 
bilities, in the form of a liquid column, that move through 
the center of the tank. Such instabilities have been shown to 
be a result of unique initial conditions and are not represen- 
tative of a typical liquid reorientation. A further discussion 
of the occurrence of these instabilities can be found in Chapter 
IV. 
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The conditions for each test are listed in Table II- 1. 

The acceleration specified is the axial acceleration for a full- 
size tank. In each case, a lateral acceleration with a magni- 
tude of approximately ten percent of the axial acceleration was 
also applied. 

Initially, the liquid was at rest, positioned as shown in 
Figure II-l. The gas/liquid interface was flat and perpendicu- 
lar to the vertical reference axis. The specified acceleration 
was continuously applied and the liquid reoriented to the oppo- 
site end of the tank. 



Figure H-l Full-Size Tank 
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Table II- 1. Test Matrix 


Test 

Number 

Liquid Volume 
(Percent) 

Full-Size Tank Acceleration 
(g) 

Tank Orientation Angle 
(Degrees) 

1 

25 

.02 

0 

2 

50 

.02 

0 

3 

75 

.02 

0 

4 

25 

.02 

30 

5 

50 

.02 

30 

6 

75 

.02 

30 

7 

25 

.02 

60 

8 

50 

.02 

60 

9 

75 

.02 

60 

10 

25 

CM 

o 

• 

90 

11 

50 

.02 

90 

12 

75 

.02 

90 

13 

25 

.04 

0 

14 

50 

.04 

0 

15 

75 

.04 

0 

16 

25 

.04 

45 

17 

50 

.04 

45 

18 

75 

.04 

45 

19 

25 

.04 

90 

20 

50 

.04 

90 

21 

75 

.04 

90 


22 


10 


04 


0 







II-4 


B. TEST SCALING 


The following variables characterize the reorientation of 
propellant within a tank: 


P = liquid density, 

= liquid viscosity, 

= liquid surface tension, 

9 = solid/liquid contact angle, 

A = acceleration acting on liquid, 

V = velocity of the gas- liquid interface, 

V = liquid volume, 
r = tank radius. 


The force (F) and pressure (P) exerted by the liquid on the 
tank during liquid motion are the measured variables. Using 
the Buckingham pi theorem, the dimensionless parameters of 
interest can be established to be: 


V 



(Froude Number) 


7T 


2 




(Bond Number) 


7T 


3 


R = 

e 


PVr 

M 


(Reynolds Number) 
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7T 


4 


7T 


5 


F 

PAv 


P 

PAr 


Measured Variables 


With regard to the motion of the liquid, the Froude number can 
be related to the Bond number and Reynolds number as follows: 


F = f (B , R ) 
r o e 


Based on numerous propellant reorientation tests, numerical 
coefficients were established so that the above relationship 
can be expressed as follows (Reference 1) : 



Considering only the relationship between the Froude number and 
Bond number, their variation is shown graphically in Figure II-2 
It can be seen that the Froude number is constant if the Bond 
number is greater than 10, This implies that surface tension 
forces are negligible, in comparison to the inertia and gravity 
forces, when Bq is greater than 10. The factor in the 
equation accounts for viscous effects as a function of the 
Reynolds number. The variation of is shown in Figure II-3. 
If the Reynolds number is greater than 50, viscous effects are 
negligible. Therefore, for any propellant reorientation which 
has a Bond number greater than 10 and a Reynolds number greater 
than 50, scaling can be based on Froude number alone. 

As will be shown later, the above requirements for B^ and 
Rg are satisfied for the propellant reorientation conditions 
being considered here, so Froude number is the scaling para- 
meter, That is 
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where the subscript "a" refers to the full-size system and the 
subscript "m" refers to the model. Therefore, 



and V~At, where t = time. Hence, 



The above equation yields the time scaling for a selected di- 
mensional scaling and the ratio of the actual to model accel- 
erations. It is independent of the liquid properties. 

The liquid properties enter into the scaling in assuring 
that the Bond number and Reynolds number are sufficiently large. 
A large liquid density helps to make both Rg and Bq large, and 
also assures that the forces due to a given volume of liquid 
will be large. Low surface tension and viscosity are also de- 
sirable . 


An evaluation of the various methods of producing the 
scaled test conditions showed that the drop tower test facility 
would best satisfy the above requirements. Martin Marietta ' s 
drop tower provides a specified low-g environment for a period 
of up to 2.1 seconds. 


The values selected for Am and 


rjjj, in establishing the 


time scaling for the test, are somewhat constrained due to the 
drop tower test facility. In order to make the time scaling 
as large as possible (ta large), A^ should be large. Large 
values of acceleration limit the test time in the drop tower 


because uue uruy capsule must be accelerated 


til ITSSpSCt to 


the drag shield; the travel distance is fixed. A value of 
0.09g is a practical upper limit for A^ and will still provide 
1.6 seconds of model test time. 



Correction Factor 



Note: Contact Angle cs 


Hattorl 


Goldsmith & Mason 
Davies & Taylor 


Bond No . = 


pAr" 


Figure II-2. Froude/Bond Number Relationship (from Reference 1) 



) ^ 1.0 10.0 

pVr 

Reynolds No. = 


Figure II- 3. Effect of Reynolds Number (from Reference I) 
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In selecting rjj^, it should be as small as possible in order 
to make tg large. However, it must be realized that the forces 
produced by the liquid within the model are proportional to the 
mass of liquid. A value for rm of 6.35 cm (2.5 inches) was se- 
lected as a suitable compromise between the two requirements. 
Therefore, the time scaling for ra = 1.5 meters (5 feet) and 
Ag = 0 .04g is 


t 

a 

t 

m 


V 0.09 60 

0.04 2.5 


= 7.35 


For a model test time of 1.6 seconds, the actual time simulated 
is 11.8 seconds. In simulating the case where Ag is 0.02g, the 
model acceleration is 0.045g and the ratio tg/tm remains the 
same. A model test period of 2.1 seconds was available at this 
lower acceleration, so ta is equal to 15.4 seconds. The desired 
accelerations were achieved with spring motors producing 133N 
(30 Ibf) and 66. 7N (15 Ibf) acting on a drop capsule weighing 
162 Kg (358 Ibm) . 


FC-43, a very dense fluorocarbon solvent, was selected as 
the test liquid. The properties of FC-43 (Reference 2) along 
with those of the actual liquid, oxygen, are listed in Table 
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TABLE II-2. LIQUID PROPERTIES 



FC-43 
at 68°F 

Oxygen at 162°R 

Density, gm/cc (Ibm/ft^) 

1.905(118.9) 

1.14(71.3) 

Surface Tension, dynes/cm 


-A 

(Ibf/ft) 

16.7(1. 14x10" -^) 

13.5(9.25x10 ^) 

Viscosity, CP (Ibm/ft-sec) 

6.5(4.36x10"-’) 

0.195(1.31x10 

Contact Angle, degrees 

0 

0 


The Bond number and Reynolds number for both the model and 
actual tank are listed in Table II-3. It can be seen that the 
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requirements of the scaling analysis have been satisfied in both 
the actual case and the model. The Reynolds number is consider- 
ably greater than 50 in both cases, so viscous effects will be 
minimal. Both Bond numbers are greater than 10; however, there 
is a large difference between the Bond numbers of the full size 
tank and the model. When the Bond number is greater than 10, 
the effect of surface tension on the motion of the liquid is 
small. The minimum model Bond number of 208 should insure this 
to be true. The lateral acceleration acting on the model tank 
moves the liquid up one side of the tank, further reducing the 
effect of surface tension on the liquid motion. The large dif- 
ference in Bond number between the full-size tank and the model 
indicates that the interface in the full-size tank would break 
up more than was observed in the model, but the general manner 
of the liquid motion will be the same. 


TABLE II- 3. BOND AND REYNOLDS NUMBERS 


Parameter 

Oxygen 

FC-43 

Bond Number 



Ag = 0.02g 

3.85x10^ 

208 

Ag = 0.04g 

7.70x10'^ 

407 

Revnolds Number 



Ag = 0.02g 

2.70x10^ 

1.73x10^ 

Ag = 0.04g 

4.13x10^ 

2,63xl0'^ 


Since complete reorientation of the liquid during the test 
was desirable, this was also considered in selecting the model 
test conditions. It was anticipated that the liquid would be 
reoriented but some oscillation of the liquid about its equili- 
brium position would still be present at the end of the test. 

Forces and pressures measured in the test may be scaled to 
the full-size tank, using the previously presented dimension- 
less parameters, as follows: 
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= 2.71 X lO""^ 

r 

a 


and 


= 0.156 

a 


C. TEST SYSTEM DESCRIPTION 

A test system that can produce the required subscale model 
test conditions and measure the liquid forces was designed and 
built for Martin Marietta's Drop Tower Test Facility. Flexi- 
bility to duplicate the varied test conditions, and sensitivity 
to record the small liquid forces, were the key requirements in 
designing the test system. 

1. Test Module - The test module consists of the tank, 
force measuring links and slider mechanism. This module is 
shown mounted on the drop capsule in Figure II-4. Figures 
II-5 and II-6 present front and back views of the box in which 
the force links and tank are mounted. 

The model tank is made of clear plastic. The domes were 
blown from sheet plastic and the barrel section was cut from a 
tube. The flange around the tank provides structural strength 
and permits the tank to be mounted at the proper angle within 
the force link yoke. Two screws in the barrel section of the 
tank allow it to be filled and drained. 

Three force links, two vertical and one lateral, allow 
all forces acting on the tank to be measured. Bearings at each 
end of the links permit only forces along the link axis to be 
measured. The bearings that are mounted on the box are self- 
aligning . 

Three flexures, perpendicular to the plane of the force 
links, prevent any motion of the tank out of that plane. The 
spring constant of these flexures is small in comparison to 


Lateral Spring Motor 
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Figure II-4. Test Module Mounted on Drop Capsule 







QUALITY 


rigure 11-5. Camera View 
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figure ii-6. Back View of Tank and Force Links 
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the spring constant of the load cells and force links. There- 
fore, the effect of these flexures on the force sensed by the 
load cells is insignificant. 

The platform of the slider is mounted to the rails with 
three linear bearings; one under the camera and two under the 
tank. A constant force spring motor provides the lateral ac- 
celeration of the slider. A spring motor with a force of 3.34N 
(0.75 Ibf) accelerated the slider at approximately 0.02g. An 
electric solenoid was used to release the slider at the begin- 
ning of the test. 

2. Dro p Test Facilit y - The complete drop capsule is 
shown in Figure II-7. Due to the rather high accelerations 
being used, evacuation of the drag shield was not necessary. 

A simple frame was mounted over the test module rather than 
sealing the drop capsule with its cylindrical cover. The spring 
motors that provide the axial acceleration of the drop capsule 
and a crush tube are mounted on the conical base. 

The total drop test system is illustrated in Figure II-8. 
The cable from the axial spring motor is extended and secured 
to the bottom of the drag shield. After releasing the drag 
shield from the top of the 23-meter (75-foot) drop tower, the 
drop capsule is simultaneously released within the drag shield. 
At the same time the solenoid is actuated releasing the slider. 
Both the lateral and axial spring motors accelerate the test 
module throughout the drop test. The drop capsule impacts the 
drag shield, with the crush tube absorbing the impact, and the 
drag shield lands in a bin of wheat at the end of the test. 

3. Ins trumenta tion - The motion of the liquid was re- 
corded with a 16-mm Mil liken DBM- 3a camera mounted on the 
slider. The film speed was 200 frames per second. Immediately 
before the drag shield was released, the camera was started and 
it was automatically stopped when the drag shield impacted the 
whea t . 

Quartz crystal load cells (Kistler Model 912) were used to 
measure the liquid forces. These load cells have a capacity of 
2220N (500 Ibf) in tension and 22200N (5000 Ibf) in compression 
providing the capability of withstanding the impact at the end 
of the test. Peak, higli frequency accelerations of up to 160g 
have been measured at impact. Due to their high degree of 
linearity, tliese load cells are fully capable of measuring the 
small forces due to tiie liquid motion. 
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Figure II-7. Complete Drop Capsule 
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The load cells were mounted in the force measuring links. 
Low noise cables were used to feed the output of the load cells 
to charge amplifiers. The charge amplifiers were located about 
half way up the drop tower to minimize the motion of the cable 
as the drag shield falls. The amplifiers were set on long time 
constant and the most sensitive scale that could be accommo- 
dated to measure the low amplitude and low frequency forces. 
Each charge amplifier input was momentarily grounded prior to 
the test, so all forces were measured with respect to zero at 
one- g . 


The output of the charge amplifiers was fed in parallel 
to both a tape recorder and a chart recorder. In order to 
filter out the vibration induced by the camera motor, a 10 Hz 
low-pass filter was used in the amplifier for the chart re- 
corder. An end-to-end calibration of the force measuring sys- 
tem was accomplished with the fixture shown in Figure II-9. 

Known weights were suspended from the hook at various positions 
with respect to the force links, the output was recorded on tape 
and then played back on the chart recorder. 

An accelerometer (Columbia Model 302-2) was mounted on the 
slider to accurately measure the axial acceleration of the drop 
capsule. The accelerometer allowed the effect of drop capsule 
drag and piston effect due to its motion relative to the drag 
shield to be measured. It was found that these effects are 
negligible. The output of the accelerometer was handled in a 
manner similar to the load cells. A low noise cable fed the 
output to a charge amplifier and its output was recorded. The 
charge amplifier was grounded prior to the test so all accel- 
erations were measured with respect to zero at one-g. 

An attempt was made to measure the pressure of the liquid 
at various points on the tank during the reorientation. A 
Kulite miniature transducer was used with signal conditioning 
and amplification configured to measure pressures over a range 
of 0 to 0.07N/cm2 (0.1 psi) . The transducer could be mounted 
in any of three positions, so that it could always be located 
near the final equilibrium position of the reoriented liquid. 
Initially, the transducer was exposed to the ullage gas and it 
became submerged in liquid during the test. It was found that 
changes in output produced by the somewhat cooler liquid con- 
tacting the transducer was of the same order of magnitude 
expected for the pressure. No usable data was obtained from 
the pressure transducer. 
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D. DATA REDUCTION 

A lateral force (F^) and two vertical forces (F 2 and F 3 ) 
were measured and recorded during the drop tests (see Figure 
II- 6 ) . The test data was manually scaled, and converted to a 
punch card data bank. The following steps were applied to each 
set of test data by a data reduction program. F 3 , F 2 and F 3 
were converted from voltages to forces by the appropriate scale 
factors. To smooth the data, somewhat, it was linearly inter- 
polated with respect to time and a moving average digital low 
pass filter was applied to the data to remove 8 Hz to 12 Hz 
noise generated by the test support structure. Figure II-IO 
depicts the shape of the filter used. The force triad was then 
transposed into the tank triad as shown in Figure II- 11. The 
following set of equations was used to perform the transposi- 
tion. 



where subscript (I) denotes the inertial triad and subscript 
(T) denotes the tank triad. The results were plotted with 
time as the ordinate. 


To facilitate comparison between the test and analytical 
results, the test data was further adjusted. As previously 
mentioned, the force guages registered "O" in Ig prior to the 
drop. The analytical model records this one-g force as a nega- 
tive force in the Zj direction. To make the analysis and test 
results compatible, the initial zero test forces were converted 
to negative Zj forces. This allows direct comparison between 
predicted and measured force time histories. 
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III. ANALYTICAL MODEL 


A two dimensional mechanical analog has been developed to 
predict forces and moments exerted by a moving liquid on its 
container walls. In particular, the mechanical analog is de- 
signed to predict the forces exerted on an orbital spacecraft, 
such as the Space Tug, due to large amplitude propellant slosh 
initiated by small accelerations induced during docking or 
other orbital maneuvers. Knowledge of these forces is neces- 
sary in the design of spacecraft control systems and docking 
mechanisms . 

The mechanical analog portrays the liquid as a point mass 
moving on a constraint surface. The constraint surface is de- 
termined by slowly rotating the tank (analytically) in a one-g 
field; the constraint surface is the locus of liquid center of 
mass (cm) locations prescribed during the rotation; assuming 
the free liquid surface is planar. This constraint surface is 
assumed to be axis- symmetric (with the tank body axis system); 
hence, one quadrant (90°) is sufficient to describe the entire 
surface. In the mechanical analog, this constraint surface is 
approximated by piecewise continuous elliptical segments. When 
the liquid cm deviates substantially from the constraint sur- 
face, the coefficients of the ellipse are updated in order to 
bring the cm back to the constraint surface. For some tanks 
and fill volumes, it may be suitable to have a single ellip- 
soid to approximate the constraint surface while others may 
require several segments to appropriately describe it. Forces 
that the liquid exerts on the container result from inertial 
reactive forces and viscous dissipative forces. 


A. EQUATIONS OF MOTION 

The two dimensional equations of motion for the coupled 
tank/fluid system are stated in canonical first order form. 
This form of the equations of motion fall within a framework 
which will accommodate an entire spacecraft even though the 
following discussions are limited to simulation of the test 
configuration. The goal is to verify the mechanical analog 
based on test results. Appendix A addresses the extension of 
the equations of motion to the coupled spacecraf t/tank/f luid 
system. 
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Figure III- 1 depicts the relationship of the tank body and 
inertial coordinate systems used in the mechanical analog. The 
constraint surface and its elliptical approximation are also 
shown , 



FIGURE III-l. COORDINATE SYSTEMS 


l J.U 


on the fluid cm can be written as follows: 





A + f 


( 1 ) 
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where: = fluid mass, 

Vp = velocity vector of the fluid cm relative to the 
inertial frame, 

A = constraint force normal to the constraint surface 
(inertial reactive force), 

f = viscous dissipative force on the fluid cm tangent 
to the constraint surface. 


The fluid cm is constrained to move in the Y'j-Z'j’ plane on the 
constraint surface, approximated by elliptical segments. Hence, 
the fluids velocity vector must be instantaneously tangent to 
the approximated surface. This is equivalent to writing: 


(Vp - V ) . V E = 0 (2) 

c a 

where: V = velocity vector of the point on the constraint 

surface (coincident with the fluid cm) relative 
to the inertial frame, 

VE = planar gradient of the constraint surface (out- 
ward normal vector), ^E a , 3e A 

ay J az ^ 


The velocity of the fluid cm may be defined as follows: 


V = V + V "e 
F a T T 


(3) 


where: = magnitude of the fluid cm velocity relative to 

the constraint surface with the tank triad as a 
reference , 

"ep = unit vector tangent to the constraint surface. 



Differentiating equation (3) yields. 






( 4 ) 


where : A 

a 


the acceleration vector applied during the test 

resolved to the tank triad, ^ ^ + A k 

a j k 


II 



FIGURE III- 2. DETERMINATION OF e^ 


From kinematics and Figure III-2, we can develop an expression 
for eT7 as follows: 
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Since the fluid cm must move in the Yx-Z'j plane, from equations 
(6) and (7), we can write 


• V 

0 = X e^) (8) 


where x e.^, determines the sign of 9. Substituting equation 

(8) into equation (5) we have an expression for 


= — ( 6n^ X e^) X e^ 


(9) 


Now substituting equation (9) into 


W 5 Vp may be written, 



+ 






( 10 ) 


Combining equations (10) and (1), we have the general equation 
of motion of the fluid cm represented in terms of the tank 
coordinate system. 


A 

a 


+ V e 
T T 



( ( ^Nq ^ ^^) X ^^) ( ^ "f f ) 


( 11 ) 


In order to solve for we can convert to a scalar equation 
by performing a vector dot product on equation (11) with "Sij. 



V e 


0 . . . 


-f .. 



X + e^ 


. since they are perpendicular, 

. since they are parallel and in opposite 
directions , 


f = magnitude of the viscous dissipative force. 


• f)(12) 
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Solving for from equation (12) yields 

2 

-f _ V r 

^ ^ ‘ I ®T ' (C X e^) x e^)f (13) 


•] 


Note: ( x e^) x = ~ej^^ .... unit inward normal vector to 


the constraint surface. 


hence, -^.^ • (( ejj^ x e^) x e^) = e^ • (-e^^) =0 they 


are 


perpendicular . 


Equation (13) can now be rewritten as follows, 


V = — - ( e • A ) 
T Mp '■ T 


(14) 


If ex is defined in terms of its components, g and h, in the 
tank triad, equation (14) can be written, (_ A , , 

1 £x " ® J ^ ' 




(15) 


From equation (15) we note that Vx is a function of the compo- 
nents of the unit tangent vector, "ex. Both of these components, 
g and h, must vary with time in order to keep ‘e-j^ tangent to the 
constraint surface as the fluid cm moves through the tank. 
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FIGURE III- 3. ORIENTATION OF UNIT TANGENT AND OUTWARD NORMAL 
VECTORS TO THE CONSTRAINT SURFACE. DEFINITION OF /J . 


From Figure III-3, it can be seen 
■e'j; is at an angle /? with respect 
III- 2, we note that ^ = 0 + 90°, 
equations (6) and (7), 


that at any instant of time, 
to the ^Yx axis. From Figure 
hence /3 = 0 . Therefore, from 



(16) 


From Figure III- 3, we can also write, 

g = cos /? 
h = sin (i 


y[7~T^ 


1 


( 17 ) 
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Equation (15) now can be written, 




“f 


COS P 


- A sin p 
k. 


(18) 


Equations (16) and (18) are the state equations for the fluid 
cm and are numerically integrated to yield the state variables 
and P . The direction of the unit tangent vector, ’c'j, must 
be initially determined in order to begin integration of equa- 
tion (18) (i.e., P Q t=0) . In the mechanical analog is 
always assumed zero. 



FldURE III-4. INITIAL DIRECTION OF 


Figure III-4 delineates, for two initial cm locations, the 
initial direction, p^, of the unit tangent vector, ’e'j. From 
kinematics we know that the velocity vector of the fluid cm 
must initially have a component in the direction of the 
D'Alembert force, F. This is equivalent to writing: 
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The D'Alembert force, F, can be shown to have the same direction 
as the negative of the applied acceleration, 



( 20 ) 


The inequality in equation (19) can only be satisfied for 
0°^ 1 1 <90°. In the mechanical analog, if V = 90°, the fluid 
will not move since it is required to move on the constraint 
surface. Determination of the value of jio requires some addi- 
tional information which is provided by the following equation. 


e^ • V E = 0 (21) 

Equation (21) states that the tangent vector must be perpendi- 
cular to the normal vector at the initial fluid cm location. 

The simultaneous solution to equations (19) and (21), recalling 
equation (17), provides the initial value of . 

1. Fluid Force Determination - The constraint force A 

(inertial reactive) is determined by performing a vector dot 

product on equation (11) with the unit outward normal vector, 

"Sn . Where 
o 


V E 





A -f V ( • e ) + _T ( ejj * ( ( ew x e ) x e_) ) 


( °Nq • ^ ) 


°N„ • £ ) 


( 22 ) 


Note: e 


-Nr 


= 0 .... they are perpendicular, 


CM -1=0 .... they are perpendicular, 

o 


(( eN„ X e^) X e.^) = -1 .... they are parallel and 

in opposite directions. 


CNq ■ X X 
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^ = -A .... assuming that the inertial reactive 

force acts opposite the outward 
norma 1 . 


Hence, equation (22) reduces to: 



(23) 


(24) 


In equation (24) A is assumed to act opposite the outward 
normal vector to the constraint surface. If, in fact, it 
should act in the same direction as the outward normal vector, 

A will have a negative sign. For example, at t=0 , when V'j=0 , 
this condition occurs. This indicates a load relief on the 
propellant tank and represents the actual phenomena that occurs, 
provided the fluid is initially oriented due to some initial ac- 
celeration gradient and is not initially in zero-g. 

The viscous dissipative force, f, is a real unknown. Its 
characteristics are not well-defined, but previous investiga- 
tions (Reference 3) indicate that significant parameters may 
be kinematic viscosity, characteristic length, and gravita- 
tional environment. The mechanical analog represents f as a 
function of velocity, Vx, and the inertial reactive force, A . 


f = M I a| + ^ I V^l ; ^ , T, > 0 (25) 

The parameters and rj are variables input to the model. 

Their values may be approximated from test data or by consid- 
eration of the fluid properties and tank construction. 

2. Additional Equations - The radius of curvature used 
in equation (16) can be determined from the elliptical surface 
equa tion . 


2 2 
E = aY + cZ'" 


1=0, Constraint Surface 


(26) 
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Solving equation (26) for Y: 



(27) 


The radius of curvature is defined as, 



, evaluated at the 
cm location (y,z) 


(28) 


Similarly, p can be defined as follows, 



, evaluated at the 
cm location (y,z) 


(29) 


In addition to the equations developed above, some other 
position variables are desirable for programming purposes. In 
particular, the fluid cm position (y,z) in the tank body sys- 
tem is needed. The initial location (@ t=0) is determined by 
a Newton Raphson iteration on fluid volume in subroutine FLUDCG, 
based on an initial acceleration field. Once integration of 
the equations of motion begins, the fluid cm location can be 
determined from the surface equation, equation (26), and the 
instantaneous taiogent vector, Tx- Recalling the definition 
of "ex and V E, we can write: 


• V E = 0 = aY cos /? + cZ sin /? 
Y - - ~ TAN j8 


(30) 


Llierc lore , 


a 


(31) 
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Combining equations (26) and (31) , we can solve for Z. 


/ 2 

I— TAN fl + c 
_\a ' 

In equation (31 ) , p is a state variable determined by integra- 
tion while a and c are the ellipsoidal surface coefficients. 
From equation (27) we can define Y in terms of Z. Therefore, 
equations (27) and (32) define the magnitudes of Y and Z. 

Their signs are initially determined by subroutine FLUDCG and 
then determined by tracking fluid cm crossings of the Y'p and 
Zx axes . 

The polar location of the fluid cm can now easily be de- 
termined as follows (Figure III-l). 


(32) 


R = 

|y^ + Z^j ^ 

(33) 

i> = 

ARCTAN 

(34) 

The 

history 

desired output of the mechanical 
of forces exerted on the tank and 

analog is the time 
supports by the 


moving fluid. Equations (24) and (25), A and f, define the 
forces acting on the fluid cm. These forces must be trans- 
formed to those acting on the tank for comparison with test. 

In Chapter II, the test configuration was identified. It 
should be noted that the force measurement system, load cells, 
measures not only the forces exerted by the moving fluid but 
also inertial forces due to the tank and support structure 
mass. To facilitate the comparison between test and analyti- 
cal results, these inertial forces have been included in the 
mechanical analog. The forces and moment (in the tank body 
system) corresponding to those measured in the tests may be 
expressed as follows 


r 


FY = f cos p + A eN, 


A. M 
J s 


(35) 
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FZ = f sin P ^ 


(36) 


MX = (f sin /3 + \ ®Nqjj) Y - (f cos p + X ®Noj) ^ (.27) 


where: = mass of tank and support structure, 

®Nq ’ ^^oj^ ~ components of the unit outward normal 


vector, 


A ^ 


3 , Summary of Equations to be Solved 


State equations; 


''i ■ % ■ ^ ■ *k ^ 

k-h 

p 


Force equations: A, = 


= [I - ( • A^)j 

f = M |a j + j\| ; p > V ^ 0 

FY = f cos |8 + A eNoj - 


FZ = f sin fi + A eNni -AM 
^ Ok k s 


(38) 

(39) 

(40) 

(41) 

(42) 

(43) 


Additional equations: P = ABS 


MX = (f sin/? + A oNok^ Y - (f cos /? + A ®Nqj) 

■ ^tovh/2' 

, \dZ>/ 


d^Y 

dz2 


Z (44) 


(45) 



Ill- 14 



( 46 ) 


(47) 


(48) 

(49) 


B. COMPUTER PROGRAM: LAMPS (LARGE AMPLITUDE SLOSH) 

The equations presented in the previous section have been 
implemented for computer solution. The program, LAMPS, has 
been written in Fortran IV compatible with the MSFC Univac 
1108. Several subroutines from the existing FORMA (Reference 
4) library at MSFC have been used in addition to those devel- 
oped under this contract. LAMPS provides time history plots 
of forces for comparison to test data in addition to detailed 
printouts of state and position variable time histories which 
track the liquid cm as it moves through the tank. 

1. General Comments on the Com p uter Simulation - Figure 
III- 5 delineates the general motion of the liquid cm through 
the tank on the constraint surface. 

The liquid travels on elliptical segments that approxi- 
mate the constraint surface. When the cm deviates from the 
constraint surface more than an allowed distance, the ellipse 
is updated in order to return the cm to the constraint surface. 
The criteria for updating is expressed as follows: 


ABS (R(cm) - R(cs)) > 


(50) 
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Inertial 

Reference 

■ » : J 


FIGURE III- 5. MOTION OF LIQUID CM IN COMPUTER SIMULATION 

(EXAGGERATED) 


where: R(cm) = distance from the tank coordinate system origin 

to the fluid cm (Figure III- 5) 

R(cs) = distance from the tank coordinate system origin 
, to the constraint surface based on the fluid 

cm's current <i> (Figure III- 5) 
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GRIT = input criteria to the program expressing the 

percent allowable deviation from the constraint 
surface based on distance R(cs) . 

If this inequality is satisfied, the ellipse approximating the 
constraint surface is updated. The update is performed using 
the current fluid cm location and a point on the constraint 
surface in the direction of fluid motion (represented by points 
labeled "u” in Figure III-5) • In matrix notation, the updated 
elliptical coefficients, a and c, are defined as follows from 
the general equation for the ellipse (equation 26) : 

y(cm) z(cm)'| ^ (1 

y(u) z(u) J 1 1 

The applied accelerations, AYI and AZI, initiate and 
maintain the motion of the fluid cm. These accelerations may 
be input to the program as constants or as functions of time. 
The only restriction to the input values is that AYI must not 
equal zero at time zero. If AYI is zero at time zero, the 
fluid cm will not move and the program will terminate execu- 
tion . 




2. Com p onents of the Simulation Program - Figure III-6 
presents a general flow chart of program LAMPS. The function 
of the subroutines used in LAMPS is detailed below. 


SURF : Defines the constraint surface table by analyti- 

cally rotating the tank in a one-g field and storing the posi- 
tion of the fluid cm. The tank is assumed axis- symme trie , 
hence, SURF only rotates the tank through 90°. The table is 
stored as values of R(cs) for given <J> values (calls FLUDCG) . 


FLUDCG : Defines the initial fluid cm location (tank body 


SybLcm^ dSSUming 


A7T = 1 


and accountinc^ foT 


tank inclination angle 9X. FLUDCG works for general cylindri- 
cal tanks with hemi-ellipsoidal domes; i.e., in the limit 
cylindrical tanks and spherical tanks. A planar free fluid 
surface is assumed and FLUDCG moves this surface around until 
the calculated fluid volume (through numerical integration) 
equals the desired fluid volume within a given tolerance. The 
free surface is always aligned parallel to the inertial plane, 
I-J (Figure III-5) . 
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igure III- 6. General Flow Chart of LAMPS 
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ECOEF : Determines the ellipse which best approximates the 

constraint surface in the region of the fluid cm. ECOEF also 
checks to see if the ellipse needs updating based on the update 
criteria and updates if required. 

ETAN : Defines the instantaneous unit tangent vector, ex, 

to the elliptical approximation of the constraint surface. 

LOCATE : Determines the location of the fluid cm, y and z, 

after the fluid begins to move. The location is determined from 
the surface equation (equation 26) and equation(30)as shown in 
equations (27) and (32) . 

YDOT : Defines the derivatives of the state variables, Vx 

and fi , for use in the integration routine. 

RKADAM : Integrates the equations of motion (defined in 

YDOT) using a fourth, order Runge Kutta (Gill modification) 
algorithm. 

IQUAD : Fortran Function which determines the quadrant 

of /? . 

OUTpx ; Prints the results of the analytic simulation at 
the time increments specified. 

TLMPLT ; On option, plots time histories if Vx, Vx, /? , /8 , 
<t> , FY, FZ, and MX. In addition, TLMPLT plots the fluid travel 

through the tank, Z vs Y. 

In addition to the above subroutines, which were written 
under this contract, the following FORMA (Reference 4) subrou- 
tines are also used: COMENT, INV5, PAGHED, PLOTl, PLOTSS, 

READ, SMEQl, START, TERPl, TERP2 , VCROSS, VDOT, WRITE and 
ZZBOMB. 


3. In p ut Format - The input format for LAMPS is defined 
below as well as the definitions of the input and output varia- 
bles. Sample input and output as well as a program listing 
are provided in Appendix B. 
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INPUT FOWMAT 1_AS 

LAS 

1000 rttAU (A6» lAt3A6) RUNN0 »UNaME LAS 

IK (hunno.lq.ahstop) STOP las 

HLAD(12Ab) TITLEl LAS 

KLA0(12Ab) TlTLt2 LAS 

HtAD(6E10,3) XL»THtTDtPCVOL*THETAA*FUEN LAS 

kLAD(bI5) NP»NTMET»MA8I EtIPRINT*NPHINT*NPLOT LAS 

KtAu(4Eio.j) vkx.chit,ueltat»eni)T las 

WLAU(4E10.3) AYI»AZI»XMUtXNUtSMASS LAS 

call COMLNT comment CAPUS« LASI card 10 ZEROS COLUMNS 1-10 LAS 

IF (AZI .EU.999, .OR.AYI .tU.999, )CAlL Rt AD ( ACCEL » NA . NC t K 1 » 3 ) 1.AS 

IF (ntable.le.o)Call head(Table»ntable»nct»ki »2) las 

GO TO 1000 las 

LAS 

LAS 

definition of input VARIABLES LAS 

las 

PUNNO = RUN NUMbER PRINTED IN PAGE HEADING, LAS 

TlTLtl = title CARO PRINTED IN PAGE HEADING, LAS 

T1TLE2 * TUlE Card PRINTED IN PAGE HEADING, LAS 

XL * LENGTH OF PROPELLANT TANK CYLINDRICAL SECTION, (L UNITS) LAS 

TR s tank RADIUS, (L UNITS) LAS 

TD = HEIGHT OF TANK DOME FROM TOP OF CYLINDRICAL SECTION, (L UNITS) LAs 

PCVOL = percentage TANK FILL ,LE, 100, LAS 

THETAX = angle IaNK IS HOTATEO ABOUT A-AXIS IN iNERtlAL TRIAD, (DEGREES) LAS 

FDEN » propellant DENSITY, (F UNI TS*S£C**2/L UN1TS**4) LAS 

NR a NUMBER OF RADIAL INTEGRATION INCREMENTS ON TANK RADIUS FOR LAS 

OElcRMlNING tank CG» SUGGEST NR= SO, LAS 

NTHET * NUMhEP OF ANGULAR INTEGRATION INCREMENTS AROUND TANK LAS 

circumference for determining tank C6« SUGGEST NTHETa 50, LAS 

NTABLE a 0, head in a TABLE DESCRIBING CONSTRAINT SURFACF (PHI VS R) , LAS 

• hit 

OF Phi from O to 90 OEG, (N,GT,0,LE,20, ) las 

IPRINT a 1, normal printout, LAS 

2* r-ULL checkout PRINTOUT, LAS 

NPRINT a lamps will print EVERY NPRINT(TH) TIME POINT, LAS 

NPLOT a Ot NO Plots will he generated, las 

a 1* GENfRaTE time history plots of vtdot» vt*betadot*beta» las 

FY«FZ»MX AND FLUID POSITION Y VS Z, LAS 

VXX a ITERATION CUlOFF (PERCENT FLUID VOLUME) FOP INITIAL CALCULAT- LAS 

ION OF fluid cg# suggest vxx= 2,0, las 

CRir 3 UPDATE CRITERIA . PERCENTAGE DEVIATION FROM R(TABLE) ALLOWED, LAS 
IF /R(ACTUAL)-< (TABLE)/, GT, (CR1T*R(TA8LE)/100,) UPDATE, LAS 

delta r a TIME INCREMENT FOR INTEGRATING THE EQUATIONS OF MOTION, (SEC) LAS 

ENUT a TIMt CUTOFF FOR PRObRAM TERMINATION, (SEC) LAS 

AYl a APPLIED Y ACCELERATION IN INERTIAL TRIAD, (L UNITS/SEC**?) LAS 

MUST NOT EQUAL 0, L*S 

a 99 Y,*PEaD in time history ACCELERATION T ABLE ,( ACCEL ) LAS 

AZI a applied Z ACCELEPATION IN INERTIAL TRIAD, (L UN I TS/SEC**2 ) LAS 

a 999 ,, read in time history ACCclERATION Table, (ACCEL) LAS 

XMU a COEF, which relates FWICTTON FORCE TO INERTIAL FORCE, (N,D,) LAS 

XNU a COEF. WHICH RELATES FRICTION FOPCE TO CENTER OF MASS LAS 

VELOCITY. (F UNiTS*SF.C/| UNITS) LAS 

SMASb a STRUCTURAL MASS ASSUMED INERT AT CENTER OF TANK TRIAD. LAS 

SMAbS IS Tank structure mass and is used in calculating las 

FORCES FOR COmPaRISuN WITH TEST OaTA.(F UNI T S*SFC*»2/L UNITS) LaS 
ACCEL a maThJX of ACCELERATION TIME MlblORItS REaD IF EIThEk AYI OR L«S 

C AZI. EQ. 999, other VALUES OF EITHER AYI OP All WILL OVERRIDE LAS 

C TABLE YaLDES, MATRIX IS AN NA A i .COLUMN 1 IS TIME (SEC), LAS 
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c 

c 

c 
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c 
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c 

c 
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column 2 IS AYl* column 3 IS AZI (L UN 1 T b/seC»*2 ) , NA.LF.20. 
AYI must not F.QUAL 0. AT TIMt= 0. 

TABLt = CONbTRAlNT SURf-ACE TABLE NTABLE X ? . NTAhLE DEFINED IN 

CALL read, column 1 IS CENTER OF MASS LUCAHON PHI (DEOPEES)* 

column ? IS cokhesponoing distance from ta>mk twiau origin 

P (L UNITS). nTaBLE.LE,20. 


definition of OUTPUT PAHAMETEHS 


time = simulation time. (SEC) 

VTDOI = FLUID CM ACCELERATION. (L UNl TS/SEC**i! ) 

VT = FLUID CM VELOCITY. (L UNiTS/SEC) 

BETAD0T= angular velocity OF VELOCITY VECTOR. (OEG«tES/SEC) 

BETA = angle the VELOCITY VECTOR MAKES RITH THE TANK TRIAD Y AXIS 
( degrees. GE. 0..LE. 3b0. ) 

XtY.Z = FLUID CM LOCATION IN TANK TRIAD. (L UNITS) 

H = radial distance from Tank triad origin to fluid cm. (l units) 

PHI « FLUID CM LOCATION AS ANGLE MEASURED FROM TANK TRIAD Y AXIS 

TO RADIAL vector R. (DEGRFES) 

AY»AZ = APPLIED ACCELERATIONS AYI AND AZI TRANSFORMED TO THE TANK 
TRIAD, (L UNITS/SEC«*2) 

ACOtCCO* COEFS. IN ELLIPTICAL SURFACE EUUATION FOR THE ELLIPTICAL 
segment representing the constraint surface, (N.O.) 

ACU*Y**2*CC0*Z**= 1.0 

RHO a RADIUS OF GYRATION OF THE ELLIPTICAL SURFACE AT X»Y,Z. 

(L UNITS) 

TANGENTS J AND K AHE COMPONENTS OF THE INSTANTANEOUS UNIT TANGENT 

VECIOR WHICH IG the DIRECTION OF THE VELOCITY VECTOR. (N.D.) 
NORMAL a J and K ARE COMPONENTS OF THE INSTANTANT ANEoUS UNIT NORMAL 
VECfOR TO THE ELLIPTICAL SEGMENT. (N.O.) 

FYfFZ » FORCES EXERTED ON TANK SUPPORTS DUE TO FLUID MOTION AND TANK 

STRUCT. Inertial forces* in tank triad. (F units) 

MX * MOMENT FXERTEU ON TANK SUPPORTS DUE TO FLUID MOTION, IN 

TANK triad. (L UNITS*F UNITS) 

KEYl = U, model free TO UPDATE ELLIPTICAL SURFACE AT WILL. 

1. last update performed until BETA ENTERS A NEW QUADRANT, 
FLUID CM IS OUTSlOF TANGENT TO CONSTRAINT SURFACE AT 

AXIS INTERCEPT* HENCE NO UPDATE IS PERFORMED, 
a 2, last update PERFORMED UNTIL BETA ENTERS A NEW QUADRANT, 
FLUID CM EXCEED CRIT CRITERIA AND IS WITHIN 20 DEGREES 

OF AN AXIS Intercept, 

a 3* LAST update PERFORMED UNTIL BETA ENTEhS A NEW QUADRANT, 
t-LUID Cm exceeds CRIT CRITERIA BUT IS WITHIN 1 UE6 OF 
AXIS INTERCEPT. HENCF NO UPDATE IS PERFORMED. 

KEY2 X 0* NO UPDATE WAS PERFORMED AT LAST TIME POInT, 

X 1. update WAS performed AT LAST TIME POINT, 

NQ8 X QUADRANT OF THE ANGLE BETA. (I *2*3* A) 

NQF = quadrant of THE ANGLE PHI. (1*2*3, 4) 


NOTES 


1) THE units of THE OUTPUT PARAMETERS ARE DEPENDANT ON THE UNITS OF 
THE INPUT parameters. EITHER METRIC OH ENGLISH UNITS MAY BE USED, 

2) L UNITS - length UNITS* IN*FT.METEHS*CM. ETC. 

F UNITS - force units, LB*K6, GRAMS* ETC, 

ALL TIME units ARE SECONDS, 

ALL ANGLE UNITS ARE DEGREES, 
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3) ZZ80MB EWHORS- III-21 

NERROR*if subroutine FLUDCG FAILED TO CONVERGE. LAS 

2* AYI* 0, at time* 0, las 

4) Subroutines called by lamps— suRFfFLUDC6.EC0EF»ETANf locate*ydot. las 

RKAOAM.IOUAD»OUTPT*TLMPlT , las 

FORMA SUBROUTINES — COMENT. INV5»PaGEHD»PLOT1 »PL0TSS»READ»SMEQ1 » LAS 

START*TERP1»T£RP2*VCR0SS*VD0T«WRITE*ZZB0MB . las 

LAS 
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IV. TEST/ANALYTICAL CORRELATION 


This chapter presents a discussion of the test results and 
provides a comparison of analytical predictions with test data. 
Appendix C presents all measured force time histories along 
with those generated by the digital simulation (Program LAMPS) 
assuming no viscous dissipative force. In addition, Appendix C 
includes a test log which delineates any problems associated 
with individual tests and a qualitative appraisal of their 
worth. 


A. OBSERVATIONS ON LIQUID MOTION 

The motion of the liquid observed in the tests adds to the 
basic understanding of the manner in which the liquid moves 
during propellant reorientation in a tank. In some of the 
early work in this area (Reference 5), the applied acceleration 
was purely axial. When the liquid interface was initially flat, 
it was found that if the Bond number was less than 10, the 
liquid reoriented along the tank walls. If the Bond number was 
greater than 10, an instability formed in the center of the in- 
terface. This instability has the form of a cylindrical column 
that travels through the center of the tank to the opposite end. 
This phenomena was studied in further detail (Reference 6) by 
tilting the tank slightly off-axis from the applied accelera- 
tion. These tests demonstrated that the instability will join 
the wall flow when the misalignment of the acceleration and 
tank is as small as one degree. When the initial interface is 
highly curved, it was found that the central instability did 
not form over a range of Bond numbers from 3 to 450 (References 
7 and 8) . 

In the tests performed for this study, the potential for 
the formation of the central instability was present since the 
Bond number was large and the initial interface was flat. Two 
factors were introduced into the tests to prevent the insta- 
bility from fully forming. First, the tank was oriented at an 
angle to the axial acceleration for most of the tests (the tank 
was oriented axially for some of the tests). Hence, an effect 
similar to that observed by Bowman (Reference 6) was expected, 
in that the lack of symmetry causes the instability to be dis- 
placed toward the tank wall. Surface tension appears to be the 
force that causes this displacement of the instability. In 
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addition, a lateral acceleration that acted perpendicular to 
the axial acceleration was always applied. In some of the tests 
the slider did not function properly and rather than applying a 
lateral acceleration, only a short duration lateral pulse was 
applied; this pulse was sufficient to produce the desired fluid 
motion; namely, the fluid adhered to the tank wall during reor- 
ients tion . 

In every test the liquid reoriented along the tank wall, 
regardless of the tank orientation and magnitude of the lateral 
acceleration. For some tests, initial formation of an insta- 
bility could be observed (Figure IV-1,3). However, it quickly 
joined the flow of liquid along the wall and disappeared. The 
leading edge of the flow adhered strongly to the tank wall, 
following the wall as it encircled the tank. Apparently, a 
small lateral acceleration occurring as the liquid first begins 
to move is all that is required to keep the liquid moving along 
the wall throughout the reorientation. 

As the liquid began to move, the fluid interface remained 
relatively flat so the motion appeared as a rotation of the 
interface about the tank center. Once the leading edge of the 
flow reached the tank dome, the fluid interface began to ac- 
quire some curvature. The liquid, in general, overshot its 
final equilibrium position, continuing on around the tank, re- 
circulating a small percentage of the liquid. Very little 
splashing of the liquid was observed, and the leading edge of 
the liquid remained attached to the tank wall. During the test 
time available, the liquid was observed to reorient, overshoot 
the equilibrium position and come to a halt. Subsequent damped 
oscillation of the liquid about its equilibrium position could 
not be observed due to test time limitations. 

A typical test is shown in Figure IV-1. This is test num- 
ber 8 in which the liquid volume was 507o, the tank was inclined 
at 60° and the smaller of the two axial accelerations was ap- 
plied to the tank. 

When the liquid volume v/as 757=, the reorientation of the 
liquid was similar to that described above, except that the 
ullage assumed the. form of a bubble and moved to the opposite 
end of the tank. The bubble followed the tank wall as it moved. 
Its surface was highly irregular due to the flow of liquid about 
the bubble. At the end of the test the bubble had become some- 
what flattened and the center of gravity of the liquid had over- 
shot the equilibrium position. 


t 
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FigureIV-1. Typical Fluid Reorientation; Test 8, 
50% Fill, 60° Tank Inclination, Ag = .045g 
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With the manner of reorientation produced in these tests, 
a geyser did not form at the reoriented position of the liquid. 
Geyser formation was categorized in Reference 7 for purely 
axial acceleration. This is another phenomena of propellant 
reorientation that is eliminated by a slight off-axis distur- 
bance to the liquid flow. 

As the reoriented liquid meets the top of the tank due to 
a purely axial acceleration, a geyser may be formed that tra- 
vels through the center of the tank. Liquid is returned to 
the bottom of the tank by this geyser. The conditions under 
which a geyser would be expected were categorized in Reference 
7. No geyser was observed in any of the tests conducted for 
this program. The lateral acceleration caused the liquid to 
flow along one side of the tank, eliminating the joining of 
the flow at the top of the tank that causes geyser formation. 


B. DISCUSSION OF TEST/ANALYTICAL CORRELATION 

In general, it is felt that the forces measured during the 
drop testing are valid. However, due to the small liquid mass 
and low acceleration levels in the test, force resolution was a 
definite problem. A study of the test results shows that much 
better force definition was obtained with the larger axial ac- 
celeration. Lateral forces were much smaller than axial forces 
due to the 1-g axial load relief at drop initiation and the 
larger axial accelerations applied. Hence, the lateral load 
cell (#1, Figure II-6) was set at a very high sensitivity in 
order to resolve these small forces. This high sensitivity 
increased the susceptability of the lateral load cell to vi- 
brations due to bearing noise, support structure vibrations, 
etc. Some of this noise was reduced or eliminated by digital 
filtering as discussed in Chapter 1. This noise is most evi- 
dent in FY for a tank inclination of 0 , and FZ for a tank 
inclination of 90°. For other tank inclinations the noise is 
masked by the relatively large axial loads. 

The moment, MX, calculated from measured forces is highly 
suspect due to the small differences between load cells. Small 
errors can completely change the sign and character of the 
measured moment. It is felt that more meaningful comparisons 
can be made between test and analytical forces (FY, FZ) than 
between test and analytical values for MX. 
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Figures IV-2 through IV-7 present comparisons of test data 
and analytical predictions (FY, FZ) for the large axial accel- 
eration («.09g) and a tank inclination of 0°. Comparisons are 
shown for liquid fill volumes of 10%, 257o, 507. and 757.. In ad- 
dition, comparisons of measured and predicted moment, MX, are 
shown in Figures IV-4 and IV-6. Figure IV-8 shows comparisons 
for 257. liquid fill and a tank inclination of 45°. 

In general, the character of the predicted and measured 
force time histories are similar. Predicted values of the 
force time histories were run assuming several values for V 
(Chapter HI, equation 25); (i was always assumed to be zero. 

A study of the comparisons indicates that the assumed form for 
the viscous dissipative force (equation 25) is not optimal for 
the simulation. As mentioned in Chapter HI, the form of this 
force is not well known. An in-depth investigation into its 
form was beyond the scope of this study, however, indications 
are that it may be non-linear. 

Observation of the presented comparisons shows that the 
predicted force time histories, trend-wise match the measured 
time histories. Although, in general, predicted forces are of 
higher magnitude than those measured. Increasing rj reduces 
the predicted magnitude but also induces time lag which de- 
grades correlation. Two possible model improvements will most 
likely improve correlation. The first is an improved form for 
the viscous dissipative force as discussed above. The second 
is the relieving of the constraint that the fluid cm must al- 
ways follow the constraint surface. In the initial moments of 
the drop, the observed instability (Section A of this chapter) 
indicates that the liquid cm actually follows a trajectory in- 
terior to the assumed constraint surface. Incorporation of 
this capability of trajectory travel and improved viscous force 
aefinition would intuitively improve timing and force level cor- 
re la tion . 

The analytic simulation appears to work best for smaller 
liquid fill volumes. This may be intuitive as per the discus- 
sion of liquid motion in Section A of this chapter. The simu- 
lation assumes that the fluid is a point mass moving on a 
constraint surface obtained by slow rotation of the tank main- 
taining a flat liquid interface. In large fill volumes (i.e., 
757o) the liquid motion was more characterized by a moving ul- 
lage bubble. This is evident in the rigid body type forces 
that were measured as shown in Figure IV-7. Liquid motion 



IV- 6 


observed for fill volumes from 107o to 507o was more representa- 
tive of that assumed in the simulation (see Figure IV- 1), 

The liquid cm location can also be used in correlating 
test and prediction data. LAMPS provides time history cm posi- 
tion data which can be compared to scaled photographic records 
of the tests. In general, liquid cm location time history cor- 
relation was good and indicated the need for some viscous dissi- 
pative force to keep the fluid from completely circulating in 
some simulation cases. See Appendix B for sample plots genera- 
ted by LAMPS. 

Overall applicability of the analytical model is considered 
to be good. It is felt that improved correlation is possible 
with further study of the nature of the dissipative force and 
improvements in the allowed liquid trajectories. The test re- 
sults are also considered good and it is felt that the measured 
forces are valid comparators for the analytical model within a 
reasonable degree of accuracy. 
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FIGURE IV-3. TEST 13, 9X = 0°, 257= FILL, Ag = .09g 
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FIGURE IV-7. TEST 15, 0X = 0°, 75 % FILL, A, 
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V. CONCLUSIONS AND RECOMMENDATIONS 


The concept, used in the analytical model, of a point mass 
moving on a constraint surface, yields results which compare 
favorably to measured test data for fill volumes up to 507o. 

Fine tuning of the analytical model may be achieved by better 
definition of the nature of viscous dissipative forces and al- 
lowing the fluid cm to move on trajectories interior to the 
constraint surface. 

1. A study should be undertaken to determine the func- 
tional characteristics of the viscous dissipative 
force as related to fluid characteristics, fill vol- 
ume and tank geometry. 

2. The mechanical analog should be modified to allow the 
fluid cm to move interior to the constraint surface, 
dependent on internal fluid bond forces and applied 
acceleration fields. 

The mechanical analog shows promise for use in the design 
of orbital control systems and the design of docking mechanisms. 

3. The mechanical analog should be expanded to three di- 
mensions and integrated into the general spacecraft 
equations for subsequent use in control system and 
loads analyses. 

The test configuration is capable of providing insight to 
the character of liquid reorientation and the forces exerted on 
spacecraft by the moving liquid. Force definition is best for 
medium fill volumes and larger applied accelerations. Small, 
short term lateral accelerations are all that is required to 
produce liquid reorientation along the tank wall and to prevent 
the development of spout instabilities. 

4. Further testing should be conducted to build the data 
bank necessary for analytical model verification. 
Testing should include various tank geometries; ogive, 
conical, etc. In addition, scale propellant manage- 
ment devices, baffels, etc., should be incorporated 
into the tests. 
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The possible application of the mechanical analog to 
reentry trajectory studies should be investigated for 
use on the shuttle external tank and other large reen- 
try bodies which may contain significant amounts of 
prope llant . 
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APPENDIX A - EXTENSION OF THE EQUATIONS OF MOTION TO THE 
GENERAL SPACECRAFT SYSTEM 


This appendix details an approach whereby the tank/liquid 
mechanical analog may be included in a spacecraft system of 
governing equations. These equations are cast into a state 
space framework that can be automated on a digital computer to 
provide the basis for both time and frequency domain analyses. 
The methodology is presented in the form of interconnected 
bodies that have certain constraints which restrict the motion 
between the bodies. The bodies could well be a spacecraft, 
tank and liquid mass combination. The ensuing discussion des- 
cribes some of the salient points relating to both the general 
form of the equations and their constraints. 

Governing Equations 

This general form of the governing system of equations has 
successfully been employed by Martin Marietta personnel to 
simulate the dynamics of interconnected spinning elastic bodies. 
The detailed derivations have been given by Mr. C. S. Bodley 
and Mr. A. C. Park.* 

A canonical first-order coupled set of equations of the 

form 


y[ = f (y^t) 


(A-l) 


is employed in the mathematical simulation. The components of 
the state vector time derivatives will be further discussed 
for the spacecraf t/tank/liquid combination. Particular atten- 
tion will be devoted to the elemental makeup of the constitu- 
ents . 

The state equations that govern the coupled vehicle/ tank/ 
liquid motion may be depicted as describing the dynamical 


* Carl S. Bodley and A. C. Park: Response of Flexible Space 

Vehicles to Docking Impact, MCR-70-2 (Vol. I), Martin 
Marietta Corporation, Denver, Colorado, March 1970. 
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motion of three bodies moving in inertial space. Constraint 
conditions are employed to both affix the tank to the vehicle 
and to define the liquid trajectory within the tank. These 
equations are stated in the following form. 



0 







(A- 2) 


The state variables of the configuration space include 
ordinary momenta, | P j- , position and attitude coordinates ,| jS 
The vector, | i® | » contains such items as Euler angles and 
inertial position coordinates. The remaining items in the 
equations will receive additional attention throughout the 
discussion. 


For a given body, k, of the system, the component ordinary, 
momenta vector, is 


P 


k 
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Further, there exists a transformation that relates the non- 
holonomic velocities, |u^ , to generalized velocities. 
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where in (A-4) the vector non-holonomic velocities |u [ con- 
tains the three projections ‘^y>“ 2 ^k angular 
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velocity vector <<* onto the body fixed axis and the three 
projections of the reference point translational velocity 
(u, V, w)k onto the body axes. The elements of (i,j = 

1,2,3) are direction cosines; the sub-matrix |'Yf(i,j = 1,2,3) 
is an orthonormal rotation transformation relating the atti- 
tude of the body fixed axis system to the inertial frame. The 
sub-matrix, , is also a rotation transformation; however, 

it is not orthonormal since it relates vector components based 
on an orthogonal basis to those of a skew basis; namely, the 
axes about which Euler rotations are measured. 


The mass matrix for body k, appears as 
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The kth component of second part of the right hand side of 
(A- 2), {njk {p} is further identified as 
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The force/ torque vector, contains the external forces and 

torques plus any stiffness and damping force that may arise 
through connections with the other bodies making up the system. 

The constraint equations (third of A-2) are written in 
terms of the non-holonomic velocities, |u J- . The coefficient 
^b] are obtained from expressions of kinematic constraint and 
these same [b] coefficients are transposed to premultiply the 
vector I Z' I , providing constraint forces and torques. 

Kinematic Coefficients 


This subsection discusses the aforementioned kinematical 
relations involving expressions of relative and absolute velo- 
cities which lead to the form of the b] coefficients. The 
discussion will focus on two adjacent interconnected bodies. 




FIGURE A-1. TWO INTERCONNECTED BODY SYSTEM 


The origins of the body reference systems are labeled m 
and n. The portions of the body where the bodies connect are 
located and labeled p and q. 


In general, for each interconnected pair of bodies, there 
will be five (5) axis systems or coordinate bases. First, 
there will be an axis system fixed to each of the points, m, 
p, q, n (see Figure A-1), The fifth axis system is a skew or 
non- orthogona 1 basis comprising direction lines or unit vectors 
about which Euler rotations are measured. The Euler rotations 
are used to describe relative attitudes between the p and q 
frames . 

At each connection joint, there will be six (6) components 
of relative velocity (three relative Euler angle rates and three 
relative translational rates that are measured along the skew 
axes) to be expressed 
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In (A-7) Rp is a (3x3) rotation transformation relating vector 
components in the p system to components in the q system. It 
transforms from p to q. The transformation is similar, and 
note that the product (RpR^ = R^) transforms from m to q. The 
matrix, relates vector components refered to orthogonal 

axes to those refered to skew axes. The matrix is a (3x3) 
skew symmetric matrix containing components of the vector 
positioning point p from m. 
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Finally, it is pointed out that certain rows from (A-7) consti- 
tute rows of |b| and other rows of (A-7) are rows of matrix |b| 
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APPENDIX B -- COMPUTER PROGRAM LISTING AND SAMPLE INPUT/OUTPUT 


The following is a complete listing of program LAMPS and all associated 
non-FORMA subroutines. Immediately following the listing are sample input 
and output for Test 16. The output is normal and not the full checkout 
option. Following the sample printout are the plots generated for this 
test case. 
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PROGRAM LAMPS(INPUT»0UTPUT»TAPE5»lNPUT»TAPE6s0UTPUTfTAPEl»FlLMPL» UAS 10 

STAPE2) LAS 13 

C»#*»**** »•***#****#***#****♦**«***»************« ***»***************'»#*#******L AS 
C LARGE amplitude SLOSH— LAMPS —ANALYTIC SIMULATION LAS 

C - •— •• LAS 

C LAS 

C LAMPS IS A PROGRAM TO SIMULATE LARGE AMPLITUDE SLOSH OF PROPELLANT LAS 

C under the application OF LOW 6 ACCELERATION FIELDS. THE SjmULATION LAS 

C Il» 2-OIMENSIONAL, the SIMULATION CONSTRAINS THE FlUID CENTER OF LAS 

C MASS TO MOVE ON A SURFACE GENERATED BY SLOWLY ROTATING THE TANK LAS 


C (ANALYTICALLY) IN 1 6 AND TRACING THE PATH THE CENTER OF MASS LAS 
C PRESCRIBES, IN THE SIMULATION THIS SURFACE IS REPRESENTED BY PIECE LAS 
C *<ISE continuous ELLIPTICAL SEGMENTS. UPDATING OF THESE SEGMENTS IS LAS 
C performed TO INSURE THAT THE CENTER OF MASS DOES NOT DEVIATE LAS 


c substantially from the constraint surface, the program outputs 

C parameters WHICH TRACE ThE FLUID CLNTER OF MaSS AS IT MOVES IN THE 


LAS 

LAS 


C tank in addition TO PRINTING TIME HISTORIES OK THE FORCES EXERTED LAS 

C BY THE FLUID ON THE TANK, LAS 


8 ********************************************************************* *********L AS 
CODED BY RL BERRY FALL 1g 74 UNDER NAS8-30690 LAi 


C**************** •***♦**•♦***♦** *********************************************L AS 
C*********************** ***************************** ************************LAS 


C martin marietta AEROSPACE LAS 

C**************** ********************************************* ***************L AS 

C INPUT FORMAT LAS 

c ————— las 

CIOOU' read (A6«I*«3A6) RUNNO.UNAME LAS 

C IK (RUNN0.EQ,4hSToP)ST0P LAS 

C READU2At>) TiTLtI LAS 

C H£A0(l2At>) TITLE2 LAS 

c HEA0(6Eiy.3) xl*tr.td,pcvol*thetax»fden las 

C KtA0(6I5) NR»NTMET,NT*8LE.IPRINT.NPRINT»NPL0T LAS 

C REAO(4E10,3) VXX.CRIT.DELTAT.ENDT LAS 

C HEAO(4E10.3) AYI«AZItXMU*XNU.SMASS LAS 

C CALL COMEnT comment CARDS. LASj CARD 10 ZEROS COLUMNS 1-10 LAS 


C IK (A2I.E0.999,,0R,AYI,EGI,999.)CAlL PFAD(ACCEL»NA.NC.K1»3) LAS 

c if(ntable.le.o)Call heao<table.ntable.nct.ki» 2 ) las 


C 60 TO 1000 LAS 

C LAS 


C LAS 

C DEFINITION OF INPUT VARIABLES LAS 


C 

C RUNNO =* RUN NUMBER PRINTED IN PAGE HEADING. 

C TITLtl « title card PRINTED IN PAGE HEADING. 

C TITLE2 a TITlE CARD PRINTED IN PAGE HEADING. 

C XL a LENGTH OF PROPELLANT TANK CYLINDRICAL SECTION, (L UNITS) 

C TR a TANK RADIUS, (L UNITS) 

C TO a height of TANK DOME FROM TOP OF CYLINDRICAL SECTION, (L UNITS) 

C PCVOL a percentage TANK FILL ,LE, 100, 

C TMETAX a ANGLE TANK IS HOTATEO ABOUT X-AXIS IN InERTIAL TRIAD, (DEGREES) 
C FOEN a PROPELLANT DENSITY, (F UNITS*SEC**2/L UNITS**4) 

C NR a number of RADIAL INTEGRATION INCREMENTS ON TANK RADIUS FOR 
C DETERMINING TANK C6» SUGGEST NRa 50, 


LAS 

LAS 

LAS 

las 

las 

LAS 

LAS 

LAS 

A #* 

U.M3 

LAS 


LAS 

LAS 


C NTHtI 
C 


NUMBER OF AN6ULAH INTEGRATION INCREMENTS AROUND TANK 

circumference for determining Tank CG. suggest NTHtla 50. 


LAS 

LAS 


C NTABLE a 0. READ IN A TABLE DESCRIBING CONSTRAINT SURFACE (PHI VS R). LAS 

C a N, LAMPS will DEFINE AN AXIS-SYMMETRIC CONSTRAINT SURFACE LAS 


C 

C 


TO BE STORED IN A TABLE aS PHI VS R WITH NTABLE VALUES LAS 

OF Phi from 0 To 90 deg, (N,GT,0,LE.20, ) LAS 


C IPRINT a I, NORMAL PRINTOUT, 

C 2. full CHECKOUT PRINTOUT, 


LAS 

las 


omcmAL PAGE J3 
OF POOR QUALITXJ 


ooonoooooooonoono onoooooooonoooooo nonn no oo on nnnnnnnnonnn nnno no 


NPHINT * 
IMPLOr s 


VXX 

M 

CRIT 

3 

deltat 

8 

ENDT 

S 

AYI 

3 


S 

AZI 

= 

XMU 

3 

XNU 

= 

SMASS 

3 


ACCEL 


TABLE » 
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LAMPS WILL PHINT EVERY NPRINT(TH) TIME POINT. 

Of NO Plots will be generateo. 

If generate time history plots of VTOOTtVTfBETADOTfBETAf 

FYfFZfMX AND FLUlO POSITION Y VS 2. 

ITERATION CUTOFF (PERCENT FLOIO VOLUME) FOR INITIAL CALCULAT- 
ION OF FLUID CGf SUGGEST VXX= ^.0. 

UPDATE CRITERIA f PERCENTAGE DEVIATION FROM R(TABLE) ALLOWED. 
IF /R(ACTUAL)-R(TABLE)/.GT. (CRlT*R(TAGLe)/i00.) UPDATE. 

TIME increment FOR INTEGRATING THE EQUATIONS OF MOTION. (SEC) 
TIME CUTOFF FOR PROGRAM TERMINATION. (SlC) 

APPLIED Y acceleration IN INERTIAL TRIAD. (L UNI TS/SEC**2 ) 
must NOT EQUAL 0, 

999.fREA0 IN TIME HISTORY ACCELERATION TABLE. (ACCEL) 
applied 2 ACCELERATION IN INERTIAL TRIAD. (L UNITS/SEC**?) 
999,fREAD IN TIME HISTORY ACCELERATION T ABLE .( ACCEL ) 

COEF. WHICH RELATES FRICTION FORCE TO INERTIAL FORCE. (N.O.) 
COEF. which relates FRICTION FORCE TO CENTER OF MASS 
VELOCITY. (F UNITS*SEC/L UNITS) 

structural MASS ASSUMED INERT AT CENTER OF TANK TRIAD. 

SMASS IS Tank structure mass and is used in calculating 

FORCES FOR COMPARISON WITH TEST DATA.IF UNI TS*SEC**2/L UNITS) 
matrix of acceleration time histories read if either AYI OR 
AZI.E0.999, OTHER VALUES OF EITHER AYI OR aZI WILL OVERRIDE 

table Values, matrix is an na x 3 fCOLUMN 1 is time (SEOf 

COLUMN 2 is AYIf COLUMN 3 IS AZt (L UNITS/SEC**2) . NA.LE.20. 
AYI must NOT EQUAL 0. AT TIM£= 0. 

CONSTRAINT SURFACE TABLE NTABLE X 2 . NTABLE DEFINED IN 
CALL READ. COLUMN I IS CENTER OF MASS LOCATION PHI (DEGREES)* 

COLUMN 2 IS C0RRESPUN0IN6 DISTANCE FROM TANK TRIAD ORIGIN 
R (L UNITS). NTABLE. LE. 20. 


DEFINITION OF OUTPUT PARAMETERS 


TIME = simulation time. (SEC) 

VTDOr = FLUID UM ACCELERATION. {L UNlTS/SEC**2) 

VT = FLUID CM VELOCITY. (L UNiTS/SEC) 

0ETaDOT= angular velocity of VELOCITY VECTOR. (DEGREES/SEC) 

BETA *i ANGLE THE VELOCITY VECTOR MAKES WITH THE TANK TRIAD Y AXIS 
( degrees. GE. 0.. LE. 3b0. ) 

XfYfZ X FLUID CM LOCATION IN TANK TRIAD, (L UNITS) 

R * radial distance FROM TANK TRIAD ORIGIN TO FLUID CM, (L UNITS) 

PHI X FLUID CM LOCATION AS ANGLE MEASURED FROM TANK TRIAD Y AXIS 

TO RADIAL VECTOR R, (DEGREES) 


AYfAZ X applied ACCELERATIONS AYI AND aZI TRANSFORMED TO THE TANK 
TRIAD. (L UNITS/SEC**2) 

ACO.CCOs COEFS. IN ELLIPTICAL SURFACE EQUATION FOR THE ELLIPTICAL 
segment REPRESENTING THE CONSTRAINT SURl'ACE. (N.D.) 
AC0*Y**2*CC0*Z»*= 1.0 

RHO X RADIUS OF GYRATION UF THE ELLIPTICAL SURFACE AT X.Y.Z, 

(L UNITS) 

TANGENTx J and K are COMPONENTS OF THE INSTANTANEOUS) UNIT TANGENT 

VECfOR WHICH IS THE DIRECTION OF THE VELOCITY VECTOR. (N.D.) 

normal = J and K are COMPONENTS oF THE I NST ANT AN I ANtOUS UNIT NORMAL 
vector TO THE ELLIPTICAL SEGMENT. (N.O.) 

FYfFZ X EORCES EXERTED ON TANK SUPPORTS DUE TO FLUID MOTION AND TANK 
STRUCT. INERTIAL FORCES. IN TANK TRIAD. (F UNITS) 

MX X moment EXERTED ON TANK SUPPORTS [)UE TO FLUID MOTION, IN 

tank triad, (L UNITS*F UNITS) 

KEYl X 0, model free to UPDATE ELLIRTICAL SURFACE AT WILL. 

i, last update performed UNTIL BETA ENTERS « NEW QUADRANT, 
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LAS 
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LAS 

LAS 
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LAS 
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LAS 
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LAS 
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LAS 
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LAS 
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LAS 

LAS 

LAS 
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LAS 

LAS 

LAS 

LAS 
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LAS 
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las 

LaS 
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KEY2 

NQB 

NUF 


NOTES 
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FLUID CM IS OUTSIDE TANGENT TO CONSTRAINT SURFACE AT 
AXIS INTEHCtPT* HENCE NO UPDATE IS PERFORMED. 

= 2» last update PEKFORMEP until beta t-NTEHS A NEW QUADRANT, 
FLUID CM EXCEED CRIT CRITERIA AND IS WIIMIN 20 DEGREES 
UF AN AXIS INTERCEPT, 

= 3. last update performed until beta enters a new quadrant, 

FLUID CM EXCEEDS CRiT CRITERIA BUT IS WITHIN 1 DEG OF 
AXU intercept, hence no update IS PERFORMED, 

» 0. NO UPDATE WAS PERFORMED aT LAST TIME POInT, 

» 1. UPDATE WAS PERFORMED AT LAST TIME POINT, 

> QUADRANT OF THE ANGLE BETA. ( 1 3. A ) 

= QUADRANT OF THE ANGLE PHI. {1*2*3»4) 


LAS 
LAS 
LAS 

las 
las 
LAS 
LAS 
LAS 
LAS 
LAS 
LAS 
LAS 
LAS 
LAS 
LAS 

LAS 
LAS 
LAS 
LAS 
LAS 
LAS 
LAS 
LAS 
LAS 
LAS 
LAS 
LAS 
LAS 
LAS 
LAS 

bil 

tJI 

«**♦*«*** ****#******S***********#«*****************tt****************«*******L AS 

common/tank/xl*tr.td.pcvol.fmass.foen.vtank.vfluid las ‘■2d 


1) The units of the output PARAMETERS ARE DEPENDANT ON THE UNITS OF 
THE INPUT parameters. EITHER METRIC OR ENGLISH UNITS MAY BE USED, 

2) L UNITS - LENGTH UNITS. IN.FT.METE«S.CM. ETC, 

F UNITS - FORCE UNITS. LBtKG.ORAMSf ETC, 

ALL TIME units ARE SECONDS. 

ALL angle units ARE DEGREES, 

3) IIBOMU ERkORS- 

NEHR0R»1» SUBROUTINE FLUOCG FAILED TO CONVERGE. 

AYI» 0. AT TImE« 0, 

A) SUBROUTINES CALLED BY LAMPS— SURF*FLU0C6.EC0EF»ETAN. LOCATE. YDOT. 

rkadam.iquad.outpt.tlmplt • 


CUMMON/STATE/BETaO.G.H»ENO( 3) .KEY1.KEY2 
C0MM0N/VECT0R/V(2) .VDT(2) .FFflT. A{ 


3) .RH0.UHK(2) »PKK(A) 


LAS 

LAS 


30 

AU 


CUMMON/TiMESS/STARTTfOELTAT.ENOT.TtTMST LAS 50 

DIMENSION TABUE<20. 2) ♦ACCELI20. 3) LAS 60 


DATA NIT»N0T/5*<>/ LAS 70 

CALL BPLT (2HNB.2HLC) LAS 75 

Kl= 20 LAS 80 

NUla 1 LAS VO 

NU2=2 LAS VS 

VVV CALL START LAS 100 

REWIND NUl LAS 110 

REWIND NU2 LAS 115 

«*«*«*«««»«««««««•*» «•«««•««»« »«•»«•*««•*««*•«**•««»««««•«*«•*••«*»» «»«L AS 120 
INPUT DATA SECTION LAS 130 

««»a»«»»»»««»»»««»«»**«**»*»««*«««*«*»»***«*««#**«*««*****«»*«**»«**«««^A^ I AO 

LAS 150 

RtAU(NlT.1000)XL.TR»TO*PCVOL.THETAX.FDEN LAS ITO 

*****RtAD control AND INTEGRATION PARAMETERS****«**********************LAS 180 

RtAO(NlTtlOlO)NR.NTHET.NTABLE.IPRINT.NPRlNT*NPLOT LAS IVO 

READ(NIT.IOOO) VXX.CRI T.DELTAT.ENOT las 200 

*****REAU applied accelerations and FRICTION FORCE COEFS***************LAS 210 

REaO(NIT.1000)VXX.CRIT.OELTAT.ENDT las 200 

REAO(NIT.IOOO) AYI.AZI»XMU*XNU.SMaSS las 220 

LAS 230 

•****PKInT input OaTA**w********#***-»**********************************LAS 2A0 

CALL PAGEHD las 250 
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«>ihIIE(NOT.2O00) 

ITE (NOl ♦2010) XL»T«*TD»PCVOL»THETAX,FDEN»NKtNTrifcT. 


mn 


VXX*NTaBLE» 

XMUtXNU 


2lPHlNTtNPRINT,NPLOT.SMASS»CRlT»OELTAT.ENOT«AYI»A2I. 

Call coMtNT 

„ (AYI,EvJ.9Q9, ,0R,A2I.Eti.999,)CALL Rt AO ( ACCEL ♦ NA » NC 

c initialization section 

£»»'»»»»»«««»«««»«««««•«««««*••«*«««« ««««««««««««««««««<»«« 

IP= 2 

IP (IPftlNT.GT.DIPa 1 
C FIND FLUID CONSTRAINT SURFACE 
IE(NTABL£.6T.<OGOTO b 
call READ(TABLE»NTA8LE«NCT»K1,2) 

DO 6 1 = DNTABLE 

table ( I « DaTABLE ( I « 1 ) *2 . * J , 14 1 5926SA/360 . 

6 continue 

GOTO 7 

b CALL SURE (NRtNTHeT«VXX»NTABLF»TABLE»Kl ♦ IP) 

C PINO INITIAL FLUID Cfi 

7 THETA=TH£TAX/57.2957T95 

IP (THETAA.NE.999,)G0 TO 10 
READ (NIT* lOOo) GY.GZ 

gXs 0,0 
GU TO 20 

10 GXa 0.0 

GY* SIN(TmETA) 

GZa COS(THETA) 

20 CALL FLU()C6(6X*GY.62»NR»NTHET*VXX*X*YtZ*IE*D 
v«HIT£(NOT*2020) 

IF(IE.NE.0)CALL ZZBOMB ( bpL AMPS. 1 ) 

XsO.O 

C FIND elliptical SURFACE WHICH APPROXIMATES CONSTRAINT SURFACE 
CALL ECOEF (XtY*Z.ACO»CCO,NTABLE*TABLE*Kl.l.CHlT) 

C FIND applied ACCEL* IN TANK TRIAD AT T» 0, 

T= 0,0 
KtY3= 0 

KtY4= 0 

IP (AYI.N£.999.)G0 TO 30 
KEY3a I 

CALL TERPl (ACCEL( It 1 ) . T»aCCEL( 1 *2) *AYI *NA. I * I »Ki . 1 > 

IP (AYI ,EU.O.O)CALL ZZB0MB(5HLAMPS.2> 

30 IF (AZI ,NE,999. )GO TO 40 
KLYAa I 

CALL TERPl ( ACCEL (l.D . T . ACCEL ( 1 ♦ 3 ) * AZI »NA. 1 ♦ 1 »Kl . 1 ) 

40 A(l)a 0,(J 

A (2) a AZI*SIN( THETA) ♦AYIACOS (THETA) 

A (3) a AZI*C0S(THETA)-AYI*SIN(THETA) 

IF (ABS(Y) .GTf.OODGOTo 15 
IF ( r.GT.O.O.ANO.AYI .GT.O.O) Yx-l ,*Y 
IF ( Y.lT.O.O.AND.AYI.LI . 0,0) Ya-l.AY 
lb IF (A8S(Z) ,GT,.001 )GOTo If 

IF (Z.GT.O.O.AND.AYI.LT.O.O) Za-l,*Z 
IF ( Z, LT, 0, 0, ANO. A Yl, GT.O.O) Z=-1.*Z 

16 continue 

C DEFINE direction OF FLUID CO VELOCITY VECTOR 
^VTaO.O 

Call etan ( aco.cco* y.z.a (2) .acd.g.m.vt) 

c*»«**pkint Tank anq fluid volumes and fluid mass*»*****» 

Call pageho 


»K1*3) 

«*«««*»»*««»«« 


*«*»««««»««««« 


C 

c 


•NITE (NOT *2030) VT ANK * VF LUI (X » FMASS 
inD i al conditions 


LAS 260 
LAS 270 

LAS 2B0 
LAS 290 

LAS 300 
*LAS 310 

LAS 320 
*LAS 330 
LAS 340 
LAS 3b0 
LAS 300 
LAS 361 
LAb 363 
LAS 364 
LAS 305 
LAS 366 
LAS 370 
LAS 3/b 
LAS 3«0 
LAS 39o 
LAS 400 
LAS 4l0 

LAS 420 
LAS 430 

LAS 440 
LAS 430 
LAS 460 
LAS 470 
LAS 4b0 
LAS 490 

LAS 49b 
LAS 500 
LAS 510 
LAS 52o 
LAS 530 
LAb 540 

LAS 5b0 
LAS 500 
LAS 5/b 
LAS 560 
LAS 505 
LAS 590 
LAS 600 
LAS 610 
LAS 620 
LAS 630 
LAb 640 
LAS 641 
LAS 642 
LAS 643 
LAS 64H 
LAS 645 
LAS 640 

LAS 64T 
LAS 650 

LAS 6bb 
LAS 6o0 

LAS 6o0 

LAS bVO 
LAS 700 

LAS 710 
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BtTA= ATANH(h»«) i as 

U (BETA.LI .0,0)8ETA = ^,»3.U15q2654*BtTA LAb 


BfcTAO* 0ETA las 

H= ( Y*«2*Z**i!) **0.5 las 

1K(ABS<Y) ,GT.1.E-?0)(>0 To 45 LAS 

iF (if.GE.O, ) PHI* 3. 141592F54/2, LAS 


iF (Z.LT.0,)PHI= -3.141592654/?, 

UU To 46 

45 PHI= ATAN2(Z*Y) 

4b IF(PHI.LT.0.U)PHI= 2.*3.14159?6544PHI 
EN0M= ( <ACO*Y)**2*(CCO*Z)**2)**O.S 
ENO(l)= 0.0 
EN0(2)= ACO*Y/£NOM 
ENO(3)= CCO*Z/ENOM 
PmU= PMI*5/. 2957795 
IF(PHD.GT.315.,0R.PHU.LT.45.)R0 TO A 7 
IE (PHD.GT.135.,AND.PrtO,LT.225.)GO TO 47 
UK* Y 

CAX* ago 

AAX= CCO 
GO TO 48 
47 Ws Z 

CAX* CCO 
AXX» ago 


LAS 

LAS 

LAS 

LAS 

LAS 

LAS 

LAS 

LAS 

LAS 

LAS 

LAS 

LAS 

LAS 

LAS 

LAS 

LAS 

LAS 

LAS 


48 DY= <-l.*CXX*W)/(AXX*( (l,-CXX*W**2)/AXX)**0.5) LAS 

0Y2=( (-1,*CXX)/(AXX*( ( I ,-CXX*8**2)/AXX)**0.5) )-( (CXX*W)*»2/ LAS 

I (AXX**2*{ (l.-CXX*w**2)/AXX)**l,5) ) LAS 

IF (AaS(0Y2) .GT, I.E-21 )GO TO 41 LAS 


HHO» 1,E*21 
GO TO 42 


LAS 

LAS 


41 HMOii ABS( ( l,»DV4*2)**l,5/OY2) 

42 CALL VOOT (ENOt AfENOA»OUM,UUM»OUM) 

F iNRa -1 .*FMASS* (ENOA- ( VT**2/PHO) ) 
FFRTa XMO*ABS(FInR) ♦XNU*A bS(VT) 

F iNRa -l.*FMASS*(ENOA-(VT**2/RHO) ) 
V(l)a VT 
V(2)a beta 

CALL YDOT 
KEYU 0 


LAS 

LAS 

LAS 

LAS 

LAS 

LAS 

LAS 

LAS 

LAS 


KEY2* 0 las 
FZ=FFRT*H*FINR*EN0(3)-A(3) *SMASS LAS 
FY»FFRT*G*FINR*EN0(2) -A (2) *SMASS LAS 
FMX* (FFHT*H*FINR*EN0(3) ) * Y- ( FFRT*G*F INR*EN0 ( 2 ) )*Z LAS 


Nub* IOUAU(PETA) 

IF ( Y.GT.O. .AND.Z.GE.U. ) nof=i 
IF ( Y.LE.O. .ANO.Z.gT.O. ) NUF*2 
1F<Y.lT.6..AND,Z.lE.O.) N0F=3 
IF (Y.GE.O..AND.Z.LT.U.) n 0F*4 
BtTO* BETA*57. 2957795 
BE TDDsVDT (2) *57.2957795 

*<HirE(NUUT» VOT(!)jYl, t!LTnn.ac;Tn,x*Y»Z 


LAS 

LAS 

LAS 

LAS 

LAS 

LAS 

LAS 


2 RH0»GfH.ENU(2) »EN0(3) .FYfFZ»FMX.KEYl.KEY2»NQH»N(JF LAS 

IF(1PRINI.NE.2)G0T0 49 LAS 


call pageho las 

"RITE (NOT»7O01 ) LAS 


N 

*c 




49 WKITE(NU2)T»V0T(1) .VT*BETUO»BFTD»PHO*Y»ZtFY#FZ.FMX LAS 

nkt=i las 

SET OP parameters FOR INTEGRATION ROUTINE— RKAOAM LAS 

NT* 0 las 

startt* 0.0 las 


730 

740 

750 

760 

761 

762 

763 

764 

7 70 
7 B 0 
790 
800 
810 
820 
825 
821 
821 
622 

823 

824 

825 

826 
027 
828 

830 

840 

841 

842 

843 

844 

850 

860 

870 

880 

6l\i 

890 

900 

910 

920 

930 

940 

960 

960 

970 

9U0 
990 
1000 
1010 
1015 
101 7 

1 /\ J t \ 

i Ub V 

1030 

1031 

1032 

1033 

1036 

1039 

10*0 

1050 

1060 
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C******«*<»****<»«*#*<HHHHt*«**<Ht»«Mt»**»*****<HMHt*«***«»**«**<»***< 

C InTEGHATION loop — SOLUTION OF EQUATIONS OF MOTION 

»««*«»««««•««««*«««•*««*< 

333 CALL RKAUAM(<JtNT) 

NKT=nkT41 

C LOCATE fluid CG IN TANK TWIaD — Y.Z 
G«C0S(V(2) ) 

M=SIN( V(2) ) 

call locate (ACO»CCO.V(E:) fYfZ) 
c calculate reqd output parameters 

via V(l) 
dtlAa V(2) 

IF (dETA,LT,0,0)UETA= 2,*3. 1 a 1S9265A*BETA 
Ks { Y**2+Z**2) **0.5 
PHI= ATAN2(Z»Y) 

IE (PHI.lt, 0,0) PHI= 2.*3,1 A1S92654+PHI 
ENOM= ( (ACO*Y) **2* (CCO*2) **2) **0,5 
EN0(2)= AC0*Y/EN0M 
EN0(3)= CCO*Z/ENOM 
PhO= PMI*5T. 2957795 
IF(PHD.GT.315..0R,PHL).LT,A5,)G0 TO 60 
IK (PHD.GT,135.,AND.PHD.LT.225.)G0 TO 60 
a= Y 

CXX= ACO 
AAXa CCO 
GU TO 61 
Z 

CXXa CCO 
AAXa ACO 

OY* (-l.*CXX*W)/(AXX*< (1,-CXX*W**2)/AXX)**0,5) 

DY2a( (-i.*CXX)/ (AXX*( ( I .«CXX*W#*2 ) /AXX ) **0 , 5 ) )•( (CXX*W)* 

1 (AXX**2*( (1,-CXX*W**2)/AXX)**1,5) ) 

1K(ABS(DY2) ,GT,1,E-21)G0 TO 62 

Rh0= 1.E-^21 
GO (0 63 

RHO= ABS( (1,*DY**2)**1.5/0Y2) 

IK(KEY3.£Q.1)CALL TEMPI ( aCCEL ( 1 » 1 )* T f ACCEL ( I »2 ) .AYI«NAt 1 
IE (KEYA.Eu.l )CALL TEMPI ( ACCEL ( I » 1 ) ♦ T ♦ ACCEL ( I ♦ 3 ) »AZI»NA*1 
A(2)x AZI*SIN (THETA) ♦AYI*COS (THETA) 

A (3) = AZI*COS(THETA)-AYI*SlN(THETA) 

Call vdot (eno,a»enoa*oum»oum*dum) 

FINH= -1.*FMASS*(EN0A^(VT**2/RM0) ) 

FFHT= XMU*aBS(FInR) ♦XNU*ABS(VT) 

FZ=FFRT*H*FINH*EN0 ( 3) -A (3) *SMASS 
FYxFFRT*G^FINP*EN0(2)-A(2)*SMASS 

FMX=(FFRT*h*FINH*FNO( 3) )*Y-{FFRT*6*FINH*EN0(2) ) *Z 
NuB= IQUAu(BETA) 

IF (Y,GT,0,,ANO.Z.GE.O.)NOF= 1 
IF (Y,LE.0.,aND.Z.gT.0.)NQF= 2 


60 


61 


•2/ 


62 

63 


»1»K1,1) 

*l»Kl,l) 


ini«k^»8*»AND,Z.LE,0.)NQF = 

IF (Y,g£,0.,AND,Z.LT.0.)NQF= 


STOKE DATA UN NUl FOR OUTPUT 
BtTD= BETa* 57, 2957795 
at TDDxVDT (2) *57,2957 T95 

'»MITE(NUl) T,VDT (1) »VT, BE T DO ,BETD» X * Y ♦ Z ♦ R ♦ PHU * A ( 2 ) *A(3) , 

. r RHO,G»H,ENC)<2) .f:no(3)»fy*fz,fmx,keyi*key2*nqb. 
IF (IPRINT,nE. 2)GOTO 5l 

KKI TE(N0T » 700 0 I T ,VDT ( 1 ) , VT, BFTDO *H£TO ♦ X » Y ♦ Z » R » PhD ♦ A ( 2 ) , 
»CCU»RH0»G,H»EN0(2) »EN0(3) .FY,FZtFMX»KEYl ♦KEY2 *NUc),N(JF ,F1 

51 WalTE(NU2) T»VDI (1) »V1 »BETOD»BFTD,PHD»Y.Z»FY»FZ.Fmx 

KEY2=0 

IF (T,GE.ENr)T)60 TO 555 


ACOtCCO* 

NOP 

A(3) , ACO. 
NR.FFrT 


LAS 

1070 

LAS 

1000 

LAS 

1090 

LAS 

1100 

LAS 

1105 

LAS 

lllo 

LAS 

nil 

LAS 

1112 

LAS 

Iizo 

LAS 

1130 

LAS 

1140 

LAS 

1160 

LAS 

1155 

LAS 

1160 

LAS 

1170 

LAS 

1100 

LAS 

1190 

LAS 

1200 

LAS 

1210 

LAS 

1210 

LAS 

1211 

LAS 

1211 

LAS 

1212 

LAS 

1213 

LAS 

1214 

LAS 

1215 

LAS 

1216 

LAS 

121 ( 

LAS 

1210 

LAS 

1220 

LAS 

1230 

LAS 

1231 

LAS 

1232 

LAS 

1233 

LAS 

1234 

LAS 

1240 

LAS 

1250 

LAS 

1260 

LAS 

1270 

LAS 

1200 

LAS 

1290 

LAS 

1300 

LAS 

1310 

LAS 

1340 

LAS 

1350 

LAS 

1360 

LAS 

13T0 

LAS 

1300 

LAS 

1390 

tsi 

Itn 

LAS 

1420 

LAS 

1425 

LAS 

142 ( 

LAS 

1430 

l-A>? 

1440 

LAS 

14*1 

LAS 

1442 

LAS 

1443 

LAS 

1445 

LAS 

1449 

LAS 

1450 
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C CHECK TO SEE IF ELLIPSE NEEDS UPDATING 
IE ( VT.LE.0.0)6OT0 70 

IF (KEY1.EvJ,0)60 to 50 
NUO* IQUAU(BETAO) 

NON* IQUAO(BETA) 

It- (NQO.EW.NQN) betao=beta 
IF (NQO.EU.NON)GO TO J33 
50 KEYls 0 

CALL ECOEF <X*Y»Z»ACOtCCO,NTABLE. table »K lt2»CHin 
C IF KtY23 I ELLIPSE UPDATED IF NOT PETU«N AND InTEGkaTE 


LAS 1460 
LAS 14/S 

LAS 14?0 
LAS 14dU 

LAS 14S0 
LAS 1500 

LAS 1510 
LAS 1520 
LAS 1530 
LAS 1540 


BETAO* BETA 
U (KEY2.Nt,l)GO TO 333 


LAS IbSO 
LAS 15b0 


CALL ETAN(AC0,CC0.Y,/,A(?) »A(3) ♦G.HtVT) 
70 IF(VT,6T.O,0)60 To 75 


LAS 1561 
LAS 15o2 


VTxO.O 

6=-l,*6 


LAS 15o3 
LAS 1564 


H*-l 

key! 


«0 


LAS 1565 
LAS 1566 


75 


?fI9iTS:«!4:j?AETA. 


2.«3.14159265A*BETA 



65 


66 


betao=beta 

V(1)=VT 

^315,, OR. PHD. lT, 45,160 TO 65 

IF(PH0.6T,I3S,, ANO.PHu.lt, 225, 160 TO 6b 
W* Y 

Cxx» ACO 
AXX= CCO 
GO TO 66 
1 

CXX» CCO 
AXXa ACO 

DY= (-1,*CXX*W)/(AXX*( (1.-CXX«W**2)/AXX)**0,5) 

DY2= ( (-1 ,*CXX) / ( AXX* ( (l.-CXX*W*«2)/AXX)»*0,5) )-( (CXX»W)*»2/ 


LAS 

LAS 

LAS 

LAS 


67 HHO* ABS( (l.*DY«*2)**l,b/DY2) 

68 ENOM* ( ( ACU*Y)*«2. (CC0*Z) **2) •*0,5 


I! 


■NO(2)a ACO*Y/ENOM 
:N0(3l* CCO*Z/ENOM 

C IFdPRINi ,NE,2)G0T0 333 

WHKE (NUT.7O021T 

WHITE (NOT»7000)bETA»G,H»ENO (2) ,EN0(3) »HH0»VT 
60 TO 333 

C OUTPUT SECTION 

Q»««« ««»»««««»«*«» 

555 CALL OUTPT (NU1.NPHINT.ENDT) 

IF (NPLOT.GT.O)CAll TUMPLT (NKT.NU2.TH) 

Go TO 999 

C***#* **♦#****■**#*««.«***«*** **^**<m 

c format statements 

i«?8 

2000 FuHmAT (///.21X.61HLAHGE AMPLITUDE SLOSH — LAMPS — AN/ 

lAL SIMULATION. /. 21X. 1H-.5X.3H . 7X . 1 H- .///.<» 1 X . i 9HI N P U T 

2 T A./.41X. 19ri********#****«***#*) 

2010 format <//»41X.3HXL».2X.F8.2./.41x.3HTRs. 2X.F8,2./.41X.3HT0*.J 
12./.38X.6MPCV0L*.3X,F7,2,/.37X.7HTHETAX=.4X.F6,2./.39X.5HFOEf 


1600 

1605 

111 } 


LAS 

1611 

LAS 

1612 

LAS 

1613 

LAS 

1614 

LAS 

1615 

LAS 

l6lo 

LAS 

1617 

LAS 

1618 

LAS 

1620 

LAS 

1630 

HI 


bSI 

1811 

LAS 

1640 

LAS 

1650 

LAS 

1660 

LAS 

1670 

LAS 

1671 

LAS 

1672 

LAS 

16 75 

LAS 

1660 

*LAS 

1690 

LAS 

1700 

*LAS 

1710 

LAS 

17e;0 

LAS 

1 726 

LAS 

1730 

•LAS 

1740 

LAS 

1750 

•LAS 

1760 

hit 

m 

CLAS 

1790 

ALAS 

1800 

LAS 

1810 

.LAS 

1820 

♦ LAS 

1830 



B-9 

I5»/»38X*6HNTHET=*2X*I5*/#40X»4HVXXa»4X*F6,2»/LAS 
3♦47X*7HNrAbL£a♦^X. I5»/»37X»7 H!PRInT»» 2X» !&#/♦ J7X»7HNPRINTs,2X» 15 « /LAS 
4i3ttXf7HNPL0T= » 1 X » 15 t / » 38X , /HSMASS* »4X»Ell,4*/ LAS 

b»39X,5HCHlT*#4X*F6.2»/»37X,7HDELTAr=»3X.F8.3f/»39x»5HENDI*,2X,F9.3LAS 
6»/»40X»4HAriSt2X»F10.4*/»40X*4HAZI*f 2X|F i 0 . 4 , / t 4 OX ♦ 5HX MU= ,4X» LAS 

7£ 1 1 • 4 » / » 40X » t>HXNU= «4XtEll.4) LAb 

2020 Format (//• 17x»31h***initial fluid c6 location***) las 

2030 FORMAT(/////tl7X»30HrANK AND FLUID CHARACTER 1 ST I rS*/ » 1 7X ♦ l5 < ?H— )« LAS 
1/* 1 7X,14HTaNK volume* #F1 1.4,2Xt8HCU UNITS»//t 16X. ibHFLUlD VOLUMELAS 

UNITS»//»18X.13hFLUID mass* fFll.4,2X*17HF -SFCLAS 
) LAS 

7000 format ( 5X*i2£10.3,/»bX,lOEiO,3»4I5t/»5X»2ElO,3) LAS 

7001 FuHMAT (//»20X,23hFULL CHECKOUT PRINT OUT*// »5X, 19HORDFR OF PRINT OLAS 

2Ur •50HT*VTDUT*VT*ttETADOT»HETA.X*Y*Z.R*PHI*A(2) » A ( 3 > ♦ ACO ♦ * 62HCC0* RLAS 
3H0*G*H*£N0(2) »EN0 (J) » F Y* fZ . MX* KEY I # KE Y2 *N0B » NUF *F TNR « F FR T » // ) LAS 

7002 Format (//♦ bX.23HUPDATE performed at T« *F9,3./*bX,4bHNEW VALUES RELAS 

ITa(HAD) *U*H.EN0(2) *Eno(3) *RH0*VT) LAS 

LAS 


END 


1840 

18S0 

18bb 

1800 

1870 

1880 

1890 

1900 

1910 

1920 

1930 

1931 

1932 

1933 

1934 
193b 
1936 
194* 
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bUbWOUTINt RAADAM(NE0»N7 ) 

RKA 

10 

c 


RKA 

20 


COMMON /riMFbS/ START!. ULLTAT,EN[)T»T»TMST 

RKA 

30 


COMMON/ VtCTOrt/Y (2) . YDT (2) tKFRT. A (3) .RHO.UR^ (2) »PKK (4) 

RKA 

40 


CUMMON/TANK/AL.TR. TD.PCVOL.FMASS.FOEN.VTANK. VFLUIO 

RKA 

bO 

c 


RKA 

60 

r<HHt*«*«**#«*#*'t»*»****<HH»**«<HMHMHHH»«<Mt****«***'»**«<‘**'t****'*»**'<HH*«HHt**'IHH>MKA 


c 

RONGE KUTFA (GILL MOO.) INTEGRATION ROUTINE 

RKA 


c 

COOED BY CARL BOOLEY » MARTIN MARIETTA CORP. 

RKA 


c 

NEO= NUMBER OF EQUATIONS 

RKA 


c 

NTx number of delta! STEPS COUNTER 

RKA 


c 

initially must be set 0 IN calling PROGRAM 

RKA 


C*»*tt*****************»***********»******»****************<»*****<n»«*tt***HKA 



Data EPSl. F.PS2 / l.E-p. l.E-2 / 

RKA 

70 

c 


RKA 

80 


IF (NT .01. 0) GO TO 10 

RKA 

'VO 


UHK(1)« 0, 

RKA 

100 


0RK(2)» 0, 

RKA 

110 


KI-SURT (0.5) 

RKA 

RKA 

158 


PRK(2)* 1.^SQrT(0,5) 

RKA 

140 


Pr<N(4)» U,5 

RKA 

IbO 


TmST» 0.0 

RKA 

160 

10 

00 120 J X l.A 

RKA 

1 70 


OIL X J 

RKA 

180 


DO 110 1-1. NEU 

RKA 

l>.tO 


L X Y0T(1)*DELTAT 

RKA 

200 


60 TO (103.l0l.101.10b) . OIL 

RKA 

210 

101 

R^«TPRK^yiL)*(Z - QRK(D) 

RKA 

RKA 

158 

103 

Sq'tSW-**' ” 

RKA 

RKA 

240 

2b0 

10b 

R = (2 - 2,*URN ( I ) ) /b. 

RKA 

260 

107 

Y(I) IS Y(l) ♦ R 

RKA 

270 

110 

QRK(I) X URKd) ♦ 3.*R - PRK(JIL)*Z 

RKA 

280 

IF (OIL .EQ. 1 .OR. OIL .EO. 3) T » T ♦ DELTAT/2. 

RKA 

2*^0 

120 

CALL YOOT 

RKA 

300 

C 


RKA 

310 

C 


RKA 

320 

300 

Nl X NT ♦ 1 

RKA 

330 


ANT X NT 

RKA 

340 


TMST X ANf*DELTAT 

RKA 

3b0 


T » START! ♦ TMST 

RKA 

360 

C 

RKA 

370 


return 

RKA 

380 


END 

RKA 

394 
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SUBHOUTINt FTAN<AA0.CC0»Y*Z»AfB.6tH» VT) 

ETAN calculates the unit tangent VFCTUH to an elliptical surface. ETAN* 20 

ETANs GJ^HK t J= UNIT VECTOR Y DIRECTION, ETAN* JO 

Ka Unit vector z direction, etan* *o 

««*« «»*«»»««*««»«*««*««»«««««»«««««««««««««*««»««««««««««««««« «ET AN* So 

Coded by rl berry sept 197a under NAse-jobyo etan* bo 

««*««««»««« «*»««*«*«««*» *»«««»«*«««««»«**»«*«'I»««*««««««»«««»«««ETAI^« 70 

«**«**««««««»**«**««*««»**««*«**»«*«*»«***«««»»«***«*««**««««*»tt»«»»«4»»E7AN* BU 

martin marietta aerospace etan* 90 

««*««««« *»««»««»««««* «»*»*««««*««««»««««*«««««««**«*««««•«»««««««««««««£ 7AN*1 00 

** subroutine arguments ETAN*U0 

ETAN*120 

AA0= elliptical surface COEF. (Y) ETAN*130 

CCO= elliptical surface COEF, (Z) etan*i^o 

T.Za current fluid CG LOCATION. ETAN*150 

A«B= APPlIFU acceleration VECTOR C0MP0NENTS»TANN TRIAD. ETAN*170 

ACCsAj^BK ETAN*1B0 

b»Ma OUTPUT COMPONENTS OF ETAN, £TAN*190 

»«*»«»*««««««*«««*««*««*«**«««««««»»»»«***»»«««**«»*»«»»«*«««»««*»«•»•«£ 7 AN*2 00 



AG=-1, 

ETAN 

12 


Xri=-1, 

ETAN 

14 


IK (G.GE.O, )XG=1, 

ETAN 

16 


IF(H,GE.0.)AH»1, 

ETAN 

lb 


IF (ABS(Y) .QT.O.OOODGO To 10 

ETAN 

20 


HHsO.O 

ETAN 

30 


G6=l,0 

ETAN 

40 


GO TO 20 

ETAN 

50 

10 

MM=(l./( ( (CCO*Z)/(AAO*YM**2,*1,0) )**0,5 

ETAN 

60 


Gbc<l,-HH**2)**0,5 

ETAN 

TO 

20 

IF (VT.6T. 0,0)60 TO 60 

ETAN 

75 


VTa 0,0 

ETAN 

7b 


AA* -l.*A 

ETAN 

60 


BB« -l,*B 

ETAN 

90 


ACKXa AAO*Y 

ETAN 

100 


ACK2a CCO*Z 

etan 

110 


IK (ACK1.GE.O..AND,ACK2.GE.O.)GO TO 30 

ETAN 

120 



im 1^8 

30 

Gb*-1 ,*GG 

ETAN 

ISO 

AO 

ACK3» GG*AA*HH*8B 

ETAN 

160 

IF 1 

FUKCE is perpendicular TO etan USE OLD ETAN FOR NEAT DELTA T. 

ETAN 

170 


IK (ACK3.EU, 0.0) return 

ETAN 

IbO 


IF(ACk3,U, 0.0)60 TO 50 

ETAN 

190 


6= GG 

ETAN 

200 


rt= rlH 

ETAN 

210 


RtTURN 

ETAN 

220 

so 

Gs •'1,*G& 

ETAN 

230 


H= -l.*HM 

ETAN 

240 


RtTURN 

ETAN 

2S0 

bU 

G=AG*GG 

ETAN 

2S2 


M=XH*MH 

ETAN 

254 


HE r URN 

ETAN 

2Sb 


End 

ETaN 

26* 


B-12 

SUBKOUTINt locate ( A AO . CCO * BE T A » Y ♦ Z ) 
COMHON/STATE/BETAO»(j*H»£NO (3) .KEYl •KfcY2 

C LOCATE calculates ThE POSITION OP THE FLUID 

KNo.iN.- jMyssfsSrWiieiisv. 0 

»«**«*««««»««*««•«««««««»««*»««» 

COOED BY HL BEHHY SEPT 1974 


CG biVEN beta FROM 


LOC 10 
LOC do 

***LOC** lO 
LOC** do 

*» 30 
40 




martin 

ft «««*««»«««» 

subroutinl arguments 


M A R I E T 


C** 

C 

C AAOa ELI 

C LCOs EL 

C BETA® ST 

C 

C YtZa OU 

C** 

CZZBOHB-NERROR 


(0 TO d*Pl RADIANS) 


It NUO=0 OH GT 4 


ZOaltO 
Y0»l,0 

lF(Y,LT.0.)YOa-l,0 
IF(2.LT.0.)Z0a-l,0 

2N»( (AA0*G**2)/( (CC0*H)**2>, 
YN« ( (I,O-CCO*ZNftft2)/AAO)**0,5 
NUO* IQUAU(BpTAO) 

NON* IQUAU(BETA) 

NUQs NUN-NUO 
IF (NDQ.EU.O)GO TO 10 


call ZZBOMH(6HLOCATE»I ) 
10 aZ= l.O*ZO 
AY» 1,0*YO 
GO TO 60 
20 AZa I,0»ZO 
A Y 3 • 1 « 0 * Y 0 
GO TO 60 
30 AZa -l,0*ZO 
AYa 1,0*YU 
Ya AY*YN 
Za AZ*ZN 
KtTURN 
ENO 


60 


■*ft**ft*ftft»**********ftft*****ft*ftft***ft«*LOC** 

SO 

under 

NAS8-30690 

LOC** 

60 

'»*««ftftftft»ftftft*ftft**«ftft»*ft«tt*«*ftft«»»«««LOC** 

To 

'ftftftftftftftftftftftftftftftftftftftftftftftftftftftftftftttttttttvftLOC** 

BO 

T A 

AEROSPACE 

LOC** 

90 

■ft*********#**ftft****ftftftft***ftft**ft****»LOC**l 00 



LOC**llO 



LOC**UO 

(Y) 


LOC**130 

(Z) 


LOC**140 

ETAN 

WITH HORIZONTAL. 

LOC**1SO 



LOC**lbO 

It tank 

TRIAOt-INPUT OLD LOCATION 

LOC**l TO 



LOC**l(J0 



LOC**190 

:tion iuuad 

LOC**200 

^****ftftftftftftft*ftft***ft«*ftft****«ft***<nnt**L0C**2l 0 



LOC 

10 



LOC 

20 



LOC 

30 



LOC 

40 

l*CCO*G 

*ft2) ) ft*0,b 

LOC 

bO 

1 


LOC 

60 



LOC 

TO 



LOC 

80 



LOC 

90 



LOC 

100 

20 


LOC 

110 

30 


LOC 

116 

10 


LDC 

120 

30 


LOC 

125 

20 


LOC 

130 

10 


LDC 

13b 

30 


LOC 

140 

20 


LDC 

1A5 

10 


LDC 

ISO 

20 


LDC 

155 

30 


LDC 

160 

10 


LDC 

16b 



LOC 

190 



LOC 

200 



LOC 

210 



LOC 

220 



LOC 

230 



LOC 

24 0 



LOC 

2bo 



LOC 

260 



LOC 

270 



LOC 

280 



LOC 

290 



LOC 

300 



LOC 

31* 



o o o 
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PONCTION lUUAD(BETA) 

C IQUAU OETERMlNEb THE QUADRANT IN WHICH BETA (0 TO 2*PI RADIANS) 

C LIES 

C*************»**«******»*»*******<nj**,»***<nnnnnm***^nj^ni**^*^inj^^^^^^^^^ 

CODED BY HL BERRY SEPT lq/4 UNDER NAS8-30690 

martin marietta aerospace 

8Da 8ETA*S7, 29577951 

IP (bo.lt,o,,and.bd.lt.-3so,)go to 10 
IP (BO,6T.O.,AND.B0.6T,360.)60 TO lO 
DO TO 20 
10 IVI2P13 BD/360, 


IQUAO 10 
IQUAO 20 
lOUAO 30 
IQUAO AO 
lUUAD SO 
lUUAD 60 
IQUAD ^0 
IQUAO «0 
IQUAO 90 
IQUAO 91 
IQUAO 92 


BO= BO-<N2PI*360.) 

20 IF (80.LT.0,)B0« 360, *80 
IF(BD,EQ*360.) IQUAO= 1 
IF (oD,GE.o.0.ANO,B0,LT,90.0) IQUAD= 1 
IF (BD.6E. 90,0, AND, BO.LT. 180,0) IQUA0= 2 
IF (dO,6E. 180.0, AND. BO. LT,270,0) IQUAD* 3 
IF <80, GE. 270.0, AND. BD.LT. 360,0) ICJUAD* 4 

return 

END 


IQUAO 95 
IQUAO 96 
IQUAO 97 
IQUAOlOO 
IQUAOIIO 
1UUAD120 
IQUA0130 
IQUAOIAU 
iQUAOlb* 



o o oo on o o on no no on ononooooooooooo 
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SUBKOUTINt ECOEF (XB»YB,ZB*ACOtCCO»NTA6LE»TABL£»Kh(» IH0D£*C«IT) ECO 10 

C0MM0N/SrATE/BETA0»G«H«EN0(3) ,KEYlfKt'Y2 ECO la 

COMMON/VtCTOH/Y (2) tYDT (2) *FFRT» A ( 3 ) ♦RHOtUMK (2 ) * PRK ( 4 ) ECO lt> 

DIMENSION TABLE<KR*1) ♦YZ(2.2) #YZI(c»2) ECO 

Q******#***«##******'**«»******'**#********-«*****»******o**«*'in»*<»***'tt*****EC0*''* 10 
C ECOEF determines THE ELLIPTICAL SURFACE COEFFICIENTS ACO AND CCO ECO** 20 

C WHICH best OtFINE the SURFACE OF ALLO ABLE FLUID Cu LOCATIONS IN ECO** 30 

C THE REGION OK THE CURRENT FLUID CG LO^-ATION- XB»YB»ZB ( lMODt»l ) .FCOFFECO** *0 
C USES THE surface TABLE GENERATED IN SUBROUTINE SURF, WILL ALSO ECO** SO 

C CHECK FLUID CG TO SEE IF SUFFICIENT DEVIATION FROM THE SURFACE ECO** 60 

C HAS OCCURED FoR ELLIPSOIDAL COEF. UPDATE TO TAKE PLACE PRIOR TO PER- ECO** 70 
C forming update (IM00E*2). 


C**UPDATE criteria- 1) R(ACTUAL) .GT, (CRIT*R(TA8LE)/100) 


C 

C 

C 

C 

C 

C 


2) CG WITHIN 10 PERCENT*R(TANG£nT) OF TANGENT 

3) R(ACTUAl) .GT.R(TANGENT) 

IF R(ACTUAL) .GE.R(TaNGENT) UPDATE PERFORMED WITH MIRROR IMAGE OF 
CURRENT CG LOCATION, IF PHI WITHIN 20 DEG OF TANGENT AXIS INTERCEPT 

OK within 10 percent INSIDE TANGENT, AXIS INTERCEPT POINT IS USED 
FOR UPDATE, 


KEYS SET BY SUBROUTINE- KEYl= N, FINAL UPDATE UNTIL FLUID BETA 

enters NEW QUADRANT 
* 0, update should continue AS REOD, 
KEY2= I, SURFACE UPDATE WAS PERFORMED, 

» 0, NO SURFACE UPDATE WAS PERFORMED, 


ECO** BO 
ECO** B1 
ECO** 82 
ECO** B3 
ECO** B4 
ECO** bb 

ECO** B6 
ECO** B7 

I68*S 98 

ECO** 91 
ECO** 92 
ECO** 9J 
ECO** 94 
ECO** 95 
ECO** 9o 

ELLIPTICAL SURFACE- ACO*X**2* AC0*Y**2*CC0*Z**2 » 1 ECO**lOO 

*************************************************** *«*****************»ECO** 110 
CODED BY HL BERRY MAY 1974 UNDER NASB-30690 £CO**120 

»*******„***#**#*^*************************»***#****#************im»**<»*ECO** 130 
•«*«»«*««»«*«*»•«••»«»•*«*«««»»«««««•««**•««•»*»*«**** **«««***»«*««««»*ECO** 1 3b 
MARTIN MARIETTA AEROSPACE EC0**l3b 

»«««**««•»««•*«««»*»«»«»»»«« »««««»» ««*«*•*«•««•**«****«•*««***»*« «««*««EC0** 13/ 
•* SUBROUTINE ARGUMENTS ECO**140 

— ECO**IbO 

ECO**lb0 
ECO**170 

ECO**1BO 
ECO**190 

ECO**200 
ECU**210 

ECO**220 
ECO**2JO 

2, compare current FLUID CG aGAInST ALLOWABLE SURF. IF OFF ECO**2<»0 
BY CRIT UPDATE ACO,CCO,.,lF NOT RETURN TO CALLING PROGRAM. ECO**2aO 
PERCENTAGE DEVIATION FROM R(TABLE) ALLOWED, ECO*»2bO 

IF H(ACTUAL) .GT, (CRIT*R(TABLE>/100) UPDATE ACO, CCO. ECO**270 


XB,YB,Z8a current FLUID CG LOCATION, 

Aco= output elliptical surface coef. 

CCOx OUTPUT elliptical SURFACE COEF, 

ntable* number of values in matrix table, 

tables MAIHIX of SURFACE LOCATIONS (/PHI/,R), SIZE <NTabLE»2) , 
KK* HOW dimension OF TABLE IN CALLING PROGRAM-, GE.NTABLE, 


DOEGs 

IMOOE= 


CKIT» 


increment to phi for Table lookup-. le. 4b degrees, 
I, determine ACO, CCO. 




ZZB0MB...NERK0R« If KR.LT.NTaBLE 




► •290 


t 1 C* C i O /\i I T V k ^ 
JWOr'V/Vi f 4 


_ TEHPc f I N V5 , ZZBOMB , WRITE 

««*•*•«»««*»•*«•«*«»*««»««««»•«««««««««*»•»««««««««»«*«**»««•»««*«« ««»«ECU**3 I U 


IF (KR.LT.NTABLE)GO to 99 

ECO 

30 

IF (KEYl.NE.O)RETURN 

ECO 

40 

DOEG= 15.0 

ECO 

bO 

Rs (XB**2^YB**2*ZB**2)**0.5 

ECO 

60 

RXY* (XB**2.Yb**2)**0,5 

ECO 

70 

IF(RXY.GT,1,£-20)GO TO 5 

ECO 

71 

PHI= 3,141592654/2. 

ECO 

/2 

GU TO 6 

ECO 

73 

PHIs ABS(ATAN(ZB/RXY) ) 

ECO 

80 

IF(IM0DE.EQ,I)G0 to 200 

ECO 

90 


o o n 


B-15 


HMD* Pril*360,/6.2831tt 
C DETEHmI t gUAOHANT OF FLUID C<j. 

IF (Y8.6T.U.,ANL).ZB.Gt;.0.)NQF= 1 
IF (YB.LE.U..AND*Z8.GT.0,)NQF= 2 
IF(YB.LT.O,.ANO.ZR.Lt.O.)NQF= 3 
IF(YB.GE.0,.AND.ZB.LT.0.>NQF* 4 
C DETEHMINE fluid DIHECTION., .quadrant of bETA. 

BtTA= Y(<d) 

NUB* IQUAD(BETA) 

C determine axis intercept angle GAFED on NOF and NDd (PrtlU) 

GO T0(10t<?0f30»40) *NUF 
10 GO T0(98»30f98»60) fNUe 
20 Go T0(50»98.60»98) tNOB 
30 6U TO(98*bO*98tSO) *NQB 
40 GO T0(60*98«b0*98) tNUB 
98 GO TO (91 f92*93f 94) »NQF 

91 IF (ABS(ZB) .LT..001.ANO,NQB.EQ.I)60 TO 50 
IF (ABS(ZB) ,LT..OOI.AND.N08,EQ,3)60 TO bO 
GO TO 99 

92 IF(A8S{YB) ,LT,. 001. AND. NOB, EQ, 4)00 TO 50 
IF (ABS(Yd) .LT,,001.AND.NOB,EQ.2)GO TO 60 
GO TO 99 

93 IF (ABS(Zd) ,LT.,001.AND.NOd,£Q,l)GO TO 60 
IF (ABS(ZB) .LT..001.AND.NQ8.EQ.3)GO TO 50 

94 iI^CaSs??!!) .LT.. 001. AND. NQB.ECJ. 4)60 TO 60 

tF<ABS(YB),LT..001.AND,NQB,EQ.2)GO TO 50 
6u TO 99 

50 PmIU* 90. 

GO TO 70 
60 PMlUa 0,0 
70 AM* -l.O 

IF (PHIU.EU.90.) AM* 1.0 
PHIUR* PHIU*6. 26318/360, 

CALL TERPB(TABLE(1*1) ♦PHIUH»TABLE ( I *2) *RU * nT ABLE » 1 ♦ I iKH » I ) 
C CHECK ON UPDATE DUE TO SURFACE DEVIATION, 

80 CALL TERP2(TA6L£(ltl) *PHItTABLE(l»2) »RP,NTABLt»l»l»KR,l) 

IF (ABS(R-WP) .LT, (CRI T*HP/100.) )GO TO 90 
IF (ABS(PHO-PHIU) ,LE,20.)KEY1*? 

GO TO 100 

C CHECK On UPDATE DUE TO R GT RTANGENT 
90 IF(H,lT, (HU/C0S(PHI) ) )GO TO 95 

C NO U^BAfE A^^UED 

95 KtY2= 0 

return 

C FLUID CG EXCEtOS CRITERIA. . .UPDATE AAO»CCO, 

100 ()NAU= DDEG*6. 28318/360. 
phin* phi ♦ ( AM*DRA0) 

IF (R.GE. (RU/COS(PhI) ) ) KEY 1=1 
IF (R.GE. (HU/COS (PHI ) ) ) KEY2*0 
iF (R.GE, (KU/COS (PHI ) ) ) return 
IF {K£Y l.£u.O)GOtO 110 

KEY1=1 OR PHI WITHIN 20 UEG OF AXIS USE PHIU FOR UPDATE. 
EXCEPT where R.GE. (HU/COS (PHin 
IN that case use CURRENT ELLIPSE 
IF (PHD.GI .89. .OH.PHD.LT. 1, )KEY1=3 
IF (PHU.6T.89. .OR.PHD.LT.l, )KEY2=0 
IF ( PHD. 6 T. 89, .OH.PHD.LT, 1 . ) RETURN 
HHIN= PHlUR 

iio Call TERP2 ( TABLE ( 1 , 1 ) .PHIN, TaplE ( 1 *2) .RN.NTABLE ♦ 1 » 1 *KP, 1 ) 
YN» RN»CUS(PHIN) 


ECO 

lUO 

ECO 

110 

ECO 

uo 

ECO 

130 

ECO 

140 

ECO 

ISO 

ECO 

IbO 

ECO 

165 

ECO 

170 

ECO 

180 

ECO 

190 

ECO 

200 

ECO 

210 

ECO 

220 

ECO 

230 

ECO 

231 

ECO 

232 

ECO 

233 

ECO 

234 

ECO 

235 

ECO 

23b 

ECO 

23 T 

ECO 

238 

ECO 

239 

m 

240 

241 

ECO 

242 

ECO 

243 

ECO 

245 

ECO 

250 

ECO 

260 

ECO 

270 

ECO 

280 

ECO 

290 

ECO 

300 

ECO 

380 

ECO 

390 

ECO 

400 

ECO 

410 

ECO 

420 

ECO 

430 

ECO 

440 

lee 

«8 

ECO 

480 

ECO 

490 

ECO 

500 

ECO 

510 

ECO 

520 

ECO 

52b 

ECO 

530 

ECO 

535 

ECU 

540 

ECO 

550 

ECO 

551 

ECO 

552 

ECO 

553 

ECO 

554 

ECO 

555 

ECO 

560 

tCO 

630 

ECO 

b'^O 
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ZNs RN*SIN(PriIN) 

ECO 

6b0 

Yi-dtU* YB**^ 

ECO 

6b0 

YZ(1.2)= ZB«*2 

ECO 

670 

YZ(2tl)= YN**2 

ECO 

6B0 

YZ(2.2)» ZN**2 

ECO 

690 

130 CALL WRlTt(YZ»2»2*2HYZ»2) 

ECO 

70U 

call lNV5(YZ»YZi,2»2) 

ECO 

710 

ACO= YZI (ItD^YZl (1*2) 

ECO 

720 

CCO* YZI (2tl)*YZl (2»2) 

ECO 

730 

CALL XRlTt(YZl»2,2»3HYZI*2) 

ECO 

7A0 

KtY2» 1 

ECO 

7b0 

hlturn 

ECO 

760 

IMODt* 1 tlNITIAL DETERMINATION OF AAO»CCO. 

ECO 

770 

200 OHAO* ODEt.*6,28318/3bO, 

ECO 

7B0 

PHIN* PHI*0RAD 

ECO 

790 

IF (PHIN.GI, (3.1A159/2.) )PH1N» PHI-ORAO 

ECO 

BOO 

GO TO 110 

ECO 

810 

99 NtRH» 2 

ECO 

820 

WHITE (6» 777) BETA t PHI. NQF*NQB 

ECO 

821 

777 F0RMAT(///*2X.2E17.4t2I5) 

ECO 

822 

IF (KR,LT.NTABLE)NERR= 1 

ECO 

830 

CALL ZZ80MB(5hEC0EF.NERR) 

ECO 

8^0 

End 

ECO 

8b* 


on on n n n n C) on 
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bUbrtOUTINt SURF(NR»NThET»VXX*NTABLE»TA8LE»KH»IP) SRF 
( UMMON/TANK/Xl* A, 0tPCVOL»FMASS»FOENfVXOXi VTOT SRF 
dimension TABLE(KR»D SRF 
NUT= 6 SRF 

C SURF COMPUTED THE LOCATION OF ThF SURFAC? ON WHICH ThE FLUID C.G, SRF** 
C IS allowed To move in a tank of DIMENSI0NS-XL»A»B with fluid fill SRF** 
C A GIVEN PERCENT (PCVOD OF TANK VOLUME. THE RE-'>ULTANT SmRFACE IS SRF** 
C stored in matrix table as 2 COLUMNS (CYLINDRICAL COORDS), SRF** 
C TABLE (NTABLE, 1) = ABS(PHI) RAOIANS SRF*» 
C table (NTABLE»2)= R LENGiH UNITS SRF** 
(;»««»« «««««»*««»«**««»««««»«««»»»«»«»««««««»«««««*«»««««« 


10 

IS 

so 

AO 

10 

20 

SO 

40 

aO 

bO 

10 

bO 


C CuOEU BY RL BERRY MAY 1974 UNDER NAS8-3069O SRF** 90 

C»*****»*«***<HHH»*******<H»**<MHt****#**-IHHt***»**************«**«*«*«#*»**Spp*«m)0 

C MARTIN MARIETTA AEROSPACE SRF**10b 

C** subroutine ARGUMENTS SRF**U0 


g 


NH = 

NTMET = 
VXX = 


NUMoEr"of“radIal integration increments on tank RADIUS, i8^**H8 

number of angular integration increments around tank CIRCUm,SRF**140 
iteration cutoff (PERCENT OF FLUID VOLUME), SRF**lb0 


NTArLEs 

table* 

kr= 

IP* 


XL* 

A = 
B* 

PCVOL* 


number of surface locations DEFINED FOR PHI* 0 TO 90 DEG. SRF**lb0 
matrix in WHICH SURFACE LOCATIONS ARE STORED-SIZE (NTABLEtZ) ,SRF**1 TO 

ROW dimension of Table in calling program, -ge.ntable, srf**ibu 

0. NO print out. SRF**190 

1 . surface location printed AND FLUDCG OUTPUT PRINTED, SRF**200 

2. surface location PRINTED. SRF**20b 

LENGTH OF TANK CYLINDRICAL SECTION, *COMMON SRF**2l0 

tank radius, *C0MM0N SRF**220 

HEIGHT OF tank DOME FROM TOP OF CYLINDRICAL SECTION. *COMMON SRF**2S0 

PERCENT TANK FILL, «COMMON SRF**240 


C 

C** 

czzbomb-nerror* i.kr.lt.ntable 

C******CALLS subroutines FLUDCG. PAGFHD.ZZBOMB 

IF (NTaBLE.LE.KR)GO to 5 
Go TO 999 
5 GA= 0, 

TA= -3.14159/2, 

DTX* (3, 14159/2, )/(NTaHLE-1) 

IPP* IP 

IF ( IP.EQ.2) IPP* 0 
DO 10 I*1.NTABLE 
04= COSOX) 

UY= SIN(TX) 

CALL FLUaCG(GX,GY,GZ.NH,NlHET.VXX,XB.Yrt.ZB. IE. IPP) 

IF (IE, NE,0)TT=TX*360,/b, 28318 

IF nE.NE,0)CALL PAGEHQ 

IE ( le.NE.O) WHITE (NOT. 1000) I.TT.6Y.6Z 

R= ( YB**2+ZB**2) **,5 

IF(A8S(YB) ,GT,i.E-20)GO TO 6 
PHI* 3.141592654/2. 


6 

7 

lu 


Gu lO 7 

PmI* ABS<ATAN(ZH/YH) ) 


[ASM 




R 


ABLE ( 
lA= ta.dfa 
If ( IP.EQ,0)PETURN 
CALL PAGEMD 

IE (NOT .1060) 
«KI 1 1 (NOT . 1 1 00 ) 




SRF* 

*2S0 

SRF* 

*260 

SRF* 

*2bb 

^SRF* 

*270 

SRF 

So 

SRF 

70 

SRF 

80 

SRF 

90 

SRF 

100 

SRF 

lio 

SRF 

120 

SRF 

130 

SRF 

140 

SRF 

ISO 

SRF 

160 

SRF 

I/O 

SRF 

17S 

SHF 

IHU 

SRF 

2U0 

SHF 

201 

SRF 

202 

SRF 

203 

SRF 

210 

SRF 

220 

SRF 

2S0 

SRF 

240 

SRF 

2S0 

SRF 

260 

bwF 

27 u 

SRF 

280 
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DO 20 I=1*NTABLE 
Tbd= TABLtn»n*360,/6.283ia 
20 »/MITE (NOr« 1200) T8R,TA0Len*2 ) 
return 

1000 FOWMAT(/»blHSUaHOUTINE SURF— FLUIOC6 FAILED TO OONVERfiE FO 
IVaLUE »13f/»5A»8HTH£TA X*.F7.S*4H 6Y**F12.S.4 h uZ=*Fl2,5) 
lOSO^FORMAT (/»20X»26HFLUIDCG CONSTRAINT SURF ACE ♦ / • 20 X » 26H 

1100 format (//» 2ox. 1 Oh/phi (deg) /»iix» ihr*/» 20 x, ioh ,io 

laoo F0RMAT(20X«F12,5,BX»F12,5) 

<?99 Call zzbomb (ahsurf* i ) 

END 


SRF 

2^0 

SRF 

300 

SRF 

310 

SRF 

320 

F SRF 

330 

SRF 

340 

— SRF 

3b0 

SRF 

3Sb 

-) SRF 

3b0 

SRF 

370 

SRF 

360 

SRF 

39* 




oo on nnnnonnn no nnnn nn nnnnn no non oo on no 
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SUBROUTINE FLUDCG (GX.GY.GZ.NR.NThET. VXX.XBAR. YBAH.ZBAR. IE. IP) 

FCG 

10 

COMMON/TANK/XL* A,R»PCV0L*FMASS*FDEN. VXOX, VTOT 

FCG 

20 

DIMENSION AX(3) 

FCG 

30 

DATA AX/lHX.lHY.lHZ/ 

FCG 

40 

real MXY.MXZ.MYZ 

FCG 

bO 

NUT* 6 

FCG 

60 




C 

C 


FLU0C6 calculates ThE LOCATION OF THE CENTEW OF UHAViTY OF FLUID FCG** 20 

IN A CYLINDRICAL TANK WITH HEM I -ELL IPSOI DAL DOMES BASED ON TH£ FCb** 30 


APPLIED acceleration VECTORS SX*oY,6Z. THE PROGHaM OPERATES BY 
USING A NEWTON/KAPHSON iteration on THE FLUID VOLUME CONTAINED 


FCG** 40 

fcg** so 


COOED BY HL BERrtY USING PROGRAM VOLUME BY J CAHhENTFP FCG** 60 

AS A model, under NAS8-3n690 MAY 197A FCG** 90 


»***««*«*»•*»«.«»«**«*««* »««««*«**»««<»«««««««««««««**»«»«*«««««»««««««« »ecG*a 100 
martin marietta AEROSPACE FCQ**li0 


** SUBHOUTINt arguments 


FCG**U0 

FCG*»140 


6X»6Y »GZ» 

NR= 

nThETs 

VXXs 

XBAKs 

YBAR* 

ZbAH* 

IE* 


IP* 


acceleration vectors applied to TANK,.Z= axial (ANY UNITS) FCG**lb0 

X»Y= lateral FCG**160 
NUMBER OF RAOIAL INTEGRATION INCREMENTS ON TANK RADIUS, FCG**170 
NUMBER OF ANGULAR INTERGRATjON INCREMENTS AROUND TANK CIRCUMFCG**160 
iteration cutoff (PERCENT OF FLUID VOLUME), FCG**190 


X COORD, OF FLUID CG LOCATION, 

YCOORO. OF FLUID CG LOCATION. 

COOkO, of fluid CG LOCATION, 

error CODE RETURNED TO CALLING PROGRAM, 

o» algorithm converged, 

1» ALGORITHM FAILED TO C0NVE«6£-CG LOCATION MAY BE WRONG, 

0. NO PRINT our, 

1* PRINT ALL FLUDCG RESULTS, 


FCG*#200 

FC6**2i0 

FCG**220 

FC6**230 

FC6**240 

FCG**2S0 

FC6*»260 

FCG**270 


** CUMMON/TAnk/ variables EXPECTFD, FCG**280 

FCG**290 

XL* length of Tank cylindrical section, fcg**300 

A* tank RADIUS, FCG**3l0 

b* height of tank dome from top of cylindrical section, FCG**320 

PCVOL* percent tank FILL.LE,100, FC6**330 

** FC6*«340 

*****CALLS SUbHOUlINE PAGEHD FCG**3&0 

** FCG**360 

*****NUTE- SUGGEST, ,, NR* bO FCG**370 


NTHET* bO 
VXX* 2.0 


FCG**3»0 

FCG**390 




M= e: FCG 60 

IF (ABS(GY) ,lT,1,E-03)M= 1 FCG 90 

IF (ABS(GZ) ,gE,1,E-03)M= 3 FCG 100 

0R= A/FLOAT(NR) FCG llO 

or* 2,*3.14159/FL0AT(NTHET) FCG 120 


9 


Nl* 20 FCG 
ICOUnT= 0 FCG 
VN= 0, FCG 
VAOXs 3.1‘»i59*A**2*XL. (4,/3.)*3,lAib9*B*A**2 FCG 
VToTs VXUX«pCVOL/l 00, FCG 
VAYs VTOr*VXX/100. FCG 
IM IP.nE.O) call PAGEnD FCG 


IF ( IP.NE.C) white (NOT» 1000)GX,OY.GZ»XL.A,8fPCVOL»NH*NTHET* VXY.AX (M)FCG 


TmAB 5= FDtN*VTuT 
Hl= A/2. 

ICOUNT* ICOUNT*! 




FCG 

FCG 

FCG 


uo 

1*0 

170 

160 

190 

191 
19* 
19b 
196 
200 

210 
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V= 0. 

MXY= 0, 

MAZ* 0« 

MYZ* 0. 

7 DO 50 

DO 50 I*i»NTMET 
K» FLOAT <J)*0R-DR/2f 
TH= FL0AT(I)*DT-0T/2. 

X= R*C0S(TH) 

Y= K*SIN(TH) 

UAs SQRT(1,-(H/A)**2) 

ZT» B*UX *XL 
ZBa -1,*<B*UX) 

XMULT* R*DR*DT 

IF (ABS(GZ) .LT.1.E-03)60 TO 493 
ZP* H1-<6X/6Z)*X-(GY/G2)*Y 
IF (ZP,LE.ZB,AND.GZ.GT,0.)GO to 60 
IF(ZP,GE.ZT.ANO.QZ.LT.O.)GO to 60 
IF (GZ.LT.O.)GO To 490 
/A= ZT 

IF(ZP.LE.ZT)Za= ZP 


V* V*(ZA-Z0)«XMUlT 
MXYa MXY*.5*IZA-ZB)*(ZA*ZB)*XMUUT 
MXZ= MXZ-»Y*(ZA-ZB)*XMULT 
MYZ= MYZ^A*(ZA-ZB)*XMULT 

GO TO 60 
490 Za« ZP 

tF (ZP.LE.ZB)Za» ZR 
V« (2T-ZA)*XMUlT 

MXY= MXY*,5*(ZT-ZA)* (ZT*ZA)*XMULT 
MXZ* MXZ*Y*(ZT-ZA)*XMULT 

MYZa MYZ^X*(ZT-2A)*XMULT 
60 TO 60 


493 IF (ABS(GY) ,lT,1.E-03)60 TO 494 
YP* H1-(GX/GY)*X 


[F (Y.GT.YP.AND.GY.GT.O, )GO TO 60 
[F(Y,LT.YP.ANO.GY,LT.O«)GO TO 60 


GO TO 10 

494 XP= H1-(GY/GX)«*Y 


IF(X.GT.AP,aND.GX.GT.O,)60 to 60 
IF <X,LT. XP.AND.GX.lt. 0.)GO TO 60 
10 V» (ZT-ZB)*XMULT 

may* MXY*.5*(ZT-ZB>*(ZT-»ZB) *XMULT 

MAZa MXZ*Y*(ZT-2 B) *XMUlT 
MY2= MYZ*X» rZT-ZB) *XMULT 

60 CONTINUE 
50 continue 

IF(ABS(V-VN).LT,VXY)G0 TO 98 
UPCa V*100./VXOX 

IF (ICOUNT.NE.1) go to 20 

) ■* U 1 
I M^l'l * — II* 

VNMla w 
Hla 

GU TO 30 

20 IF UCOUNT.6T.2)GO TO 27 
HNa HI 

Yns V 
.0 TO 30 

27 VNMla VN 
VNs V 
HNMla HN 
HNa HI 


FCO 

220 

FCG 

2J0 

FCG 

240 

FCG 

2!»0 

FCG 

250 

FCG 

270 

FCG 

250 

FCG 

290 

F CG 

300 

FCG 

310 

FCG 

320 

F CG 

330 

FCG 

340 

FCG 

350 

FCG 

350 

FCG 

3T0 

FCG 

350 

FCG 

36s 

FCG 

369 

FCG 

390 

FCG 

400 

FCG 

410 

FCG 

420 

FCG 

430 

FCG 

440 

FCG 

450 

FCG 

451 

FCG 

452 

FCG 

453 

FCG 

454 

FCG 

455 

FCG 

456 

FCG 

457 

FCG 

450 

FCG 

4I0 

FCG 

4ti0 

FCG 

490 

FCG 

495 

FCG 

500 

FCG 

510 

FCG 

520 

FCG 

525 

FCG 

530 

FCG 

540 

FCG 

550 

FCG 

560 

FCG 

570 

FCG 

58 V 

FCG 

590 

FCG 

600 

FCG 

6 1 0 

FCG 

620 

FCG 

630 

FCG 

640 

FCG 

650 

FCG 

650 

F CG 

6T0 

FCG 

680 

FCG 

690 

FCG 

700 

FCG 

710 

FCG 

720 
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30 If- ( IP,EO»0)60 TO 40 

*»PITE(N01 *1100) Hl#V*UPC*VN»VNMl*ItOUNT 
40 IF (1COUN1 .Eo.neo to <? 

HNP1= HN-(V-VT0T)»( (HN-HnmI )/(VN-VNM 1) ) 
If (ICOUNT.tiT.NDGO TO 2b 
HI= HNPl 
bu 10 9 

98 UPC= V*100,/VXOX 


FCO 730 
FC6 740 
FCG 7b0 
Fee 7fa0 
Fee 770 
Fee 780 

Fee 790 
Fee 800 


IF ( IP.NF.O) WHITE (NOT* 1200) Hi* V»UPe» I COUNT 
IF (ieoUNT,GT.Nl.ANO.IP.NE.O)WRITE<NOT* 2000 )NI 


Fee 810 
FCG 820 


AoAH» MYZ/V 
TbAHs MXZ/V 


Fee 

Fee 


ZbAHs (MXY/v)-(XL/ 2. ) Fee 
IF (IP.NE.O) WRI TE (NOT*1300) X8AR*YBAR*ZBAR FCG 
IF (IP,NE.0 )wHITE(NOT* 1400)VXOX*VTOT Fee 
it= 0 Fee 
IF (IC0UNT.GT.NI) IE« 1 FCG 

return Fee 


C*»*F0HHAT statements*** fcg 

1000 Format (//////* 50x* 18HSUBR0UTINE FLeiUCe*///* lox* IOHINPUT OATA*/*10FCe 

IX* I OH */*10X*16hTANK ORIENT AT I0N*8X * 3HGX» *F 1 2 ,5 * / , 34 X . 3HGFCG 

2Y=*F 12,5*/ *34X*3hgZ=*f 1^.5*/* 1 OX* 13hTaNK GEUMETRY* 1 1X*3MXL = *F12,5»FC6 


3Z*35X*2HA=»Fl2,5./*35X»2hb=,F12.5*/* 10X*12H» FLUID F IlL*9X,6hPCV 0LFCG 
4=*F 12,5*/* lOX* 16HINTEGRAT10N PARAMS * 6X * 3HNR= ♦ *X ♦ I3*/*3lX*6HNTHET«*FCG 


830 

040 

8bO 

880 

86b 

8/0 

88U 

090 

900 

910 

920 

930 

940 

9b0 


5*?iI 3*/*10X*48HITERATION CUTOFF 
6F12, §*/////, bX*Al *bH “ 


• w.j wy.wr. ABS(V-VN) ,LT. (VXX*VFlUIU/100.)=* 
AXIS*12X*5HFLUI0*/*7X* 9H1NTERCEPT*9X.Auv 






7E*8X*6H$ FILL*10X*4HV(N) * 12X*6HV (N-1 ) * 8X * 9HI TEHA T ION * / * 2 ( 5X * 1 2H FCG 

*5X*7H *2(5X*12H ) *5X*9H FCG 

1100 FURmaT(2(5X*F12,S) *5X*F7,3*2(5X*Fl2,5) *8X*I3) FCG 

1200 FURMAT(2(6X«F12.S) *5X*F7,3*42X*I3) FCG 

1300 FURMAT(10X«22hFLUID C6 IS LOCATED AT , / * I OX *5HX8AR»*F 12,5 */ * 1 OX * FCG 

lbMYBARs*Fi2,b*/*10X*bMZBAR«*F12,5) FCG 

1400 FORMAT (/*11X*21H*****V0LUMES - VTANK**F 12.5* /* 25X * 7HVFLUID® *F1 2,5 ) FCG 
2000 FORMAT (/*9X,20HNO CONVERGENCE AFTEH*I3»23H ITERATIONS* BEST GUES5)FCG 
ENU fcg 


960 

970 


980 

990 

1000 

loio 

1020 

1030 

1035 

1040 

105* 



SUBROUTINt Yl)OT YDOT 
COMMON/TANK/xL»OtB»PCVOL.FMASS»FDEN»VXOXtVTOT YDOT 
COMMON/VECTOR/Y (2) tYOT (^) »FFRTt A (3) tRHOtURK (2) .PRK (4) YDOT 
COMMON/STaTE/BETaO*G.H»ENO(3) ,KEY1*KEY2 YDOT 
dimension ET(3) »AI (3) YDOT 


C DEFINES derivatives FOR USE IN RKADAM YDOT 


Q****************«*«*»***#<MHMHHHt«*«******»******««HMMHHt«MHHMHHHHHHj^,<m<j<jyyQ] 


HETA=Y(2) YDOT 
Y0T(1)= -1.*FFHT/FMASS-(A(2)*C0S(BETA) )-(A(3)*SIN(BETA) ) YDOT 
ET(1)= 0. YDOT 
ET (2)= G YDOT 
ET (3)= H YDOT 
Call VCROSS(ENO«ETtAl«UUMtOUH«AIMAAtOUM) YDOT 
YOT(2)= Y(1)/RH0 YDOT 


IF(AI (1) ,LT,0.) YDT(2)= -1*YDT(2) 

return 

END 


YDOT 

YDOT 

YDOT 


10 

Is 

20 

2b 

jO 

40 

bO 

60 

65 

TO 

60 

S#0 

100 

110 

120 

130 

140 

lb* 


o o 
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SUBHOUTINt OUTPT (NUl *NPH1NT»ENDT) OUT 

DIMENSION A(E2»13) ♦IA(4,13) OUT 

OUTPT prints the OUTPUT FOR PROGRAM LAMPS OUT 

OUT 


C COOED BY HL BERRY SpPT 1974 UNDER NAS8-J0t>90 OUT 

c martin marietta aerospace out 

0 »»««»»«»«»«»»»«««««»•»«»»«»«»«««»«««««««»»«»»«»«*«*«»»«**»««««*«»«««««« Q(JT 

RE»<IND NUl OUT 

Nur= 6 OUT 

Call pagehd out 

WKI IE (NOT.^OOO) OUT 

NP= 0 OUT 

M= 1 OUT 

HEAD (NUl) (A(I»M) «I=1»22) .(IA(J*M),J=1*4) OUT 

100 DU 10 1= ItNPRiNT OUT 


10 

19 

20 


30 


40 


50 


NP= NP*1 

IF (NP.EO.NPRInT ) GO TO 20 

REA0(nU1)T,0UM 
if(t.ge.endt)go to 19 
continue 

CALL 2ZBOmB( 5M0UTPT« U 
BACKSPACE NUl 

M3 M4 1 

READ (NUl )(A(l*M)*Islf22).(IA(J*M)*J>l»4) 

IF (A( l.M) .GE.FNDT)GO TO 30 

IF(M,GE.13)G0 to 30 

NPs 0 

GO TO 100 

DO 40 KaltM 

WRITE (NOT»20l0) (A(I»K) flal* 10) 

continue 

WHI IE(NOT.2020) 

UO 50 K«i»M 

write (NOT 1 2030) A (l»K)»(A(l»K),I = U*i9) 


8 Sf } 


out 

our 

out 

OUT 

out 

OUT 

OUT 

OUT 

OUT 

OUT 

OUT 

OUT 

OUT 

OUT 

OUT 

OUT 

OUT 


T»20A0) 

DO 60 K^ltH 

WKI TE (NOTt2050) A (1»K)»(A(I*K)» 1=20*22) »(IA(J»K)*J=1*4) 

60 continue 

IF ( A ( 1 f M) .GE.E nOT) GO TO 200 

Call pagehd 

WRITE(N0T *2000) 

NPs 0 
M= 0 

GO TO 100 
200 CALL pagehd 

<*RITE (NOT *3000) 

return 

FORMATS — OUTPUT 


10 

20 

30 

40 

50 

tiO 

70 

60 

90 

100 

110 

115 

120 

130 

140 
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160 

170 


60 

90 


200 

210 

220 

230 

240 

250 

260 

2T0 

280 

290 

300 

310 

320 

330 

340 

350 

360 


88T 


OUT 

OUT 

OUT 

OUT 

OUT 

OUT 

OUT 

OUT 

OUT 

OUT 

OUT 

OUT 

OUT 


2000 FuRMaT (///♦4lX*65HA NALYTICAL RESULT S » // ♦ 1 X ♦ ¥>0H***OUT 
l<n»««i>*«nnnnni«*<nnnnMH»*#**#<M>**<MHnH»*«***ft»****<»,j2HSYSTEM STATE *54 OUT 

3t»6x,bHVTDOTf8X*2HVT * 1 OX t 7H6ET AOOT » 5X , 4HBET A * 1 0 X 1 1 HX * 1 1 X ♦ 1 HY * 1 1 X» lOUT 

4h/» 11X*1HH»10X*3HPH1*/*1X»11 (lOH ) *5H- -*/) OUT 

2010 Format (lXtF9.3*3(3X»El 0.3) »?X.F7,2*3X.4(2X»E10.3) ,2X»F7.2) OUT 

2020 format (/* 13X*22H****APPL IED ACCEl****** 3X» J3H<»********<»***SURFACF*0UT 
l**<M>««6»<nnn,,4X»lf?H»**«»*T ANGFNT*****. 6X , I 8H*****»*N0RMAL***«* ♦ / , OUT 
24A,4HTIME*9X*2HAY. 10X.2HAZ. 10X.3HACO»9X»3HCCO»9X.3 HRhu» 9X, IHj, 1 1 X*OUT 

ORIGINAL PAGE IS 
OF POOR QUALITY 


390 

400 

410 

420 

430 

440 

450 

460 

4) 0 
46(J 
490 
500 
510 
520 
530 
5‘*0 

550 

560 

5) 0 
560 
59(j 
600 
610 
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31HK#llX»lHJ*llX»lHK»/f 1X,11 (lOH »6H */) OUT 

2030 FO«MAT(lAfF9.3»lX*5(2X»ElO,3) ,4X*4(F6.3»6X) ) OUT 

2040 FOHMAT (/♦ 13X»24H****FLUID FORCES ON TANK»4 (4H***«) ♦ l4X» 1BH«*******0UT 
1KEYS******«12X»11H*QUA0KANTS*,/*4X»4HTIME»9X»6HFY *8X»6HFZ »OUT 
2sXt8HMX tl7X.4HKEYl*10X,4HKEY2» 12X#3HNQB»5X»3HNQF./» IX .27 (4H-0Ur 

3 — )«/) our 

2050 FORMAT ( 1X.F9.3, lx *3 (2X*E 12.4) . 13X. I5»9Xt I5f 1 IXt lb.3X, 15) OUT 

3000 FORMAT (//’/////»38X*8(4H*«*») tlH*»/*38X»33H**LAMPS SUCCESSFULLY TFROUT 
lMlNATED**»/»38Xt8 (4H****) ♦ 1M*) OUT 

END OUT 


620 

630 

640 

6b0 

6b0 

670 

68u 

690 

700 
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SUBROUTINE TLMPLT ( NKT »NPL » TR ) PLOT 10 

DIMENSION PLO(200»UJ »PT (4) PLOT 20 

DATA PT(1) »PT(2) tPT(3) ♦PT(4)/10HLARGE AMPL • 1 OHITUDE SLOStlOHH SIMU LOT JO 

$LATI*2H0N/ plot 40 

C PLOTS OUTPUT PLOT 46 


XS=-1,*TR 

PLOT 50 

XD=(2,*TH)/10. 

plot 60 

XNAME1S3HSEC 

PLOT 70 


PLOT 60 
plot 90 

YNAM2=3HVEL 

YNAM3»6HBETA0T 

i?S 

YNAM4a4HB£TA 

PLOT 120 

YnAMS=3HPHI 

PLOT 130 

ynam^^Ihfy 

m 15S 

YNAM8=2HF2 

YNAM9»2HMX 


rewind npl 

PLOT 180 

N£NO=200 

PLOT 190 

IF ( NKT.lt. 200) NENOsNKT 

PLOT 200 

10EL«NKT/205*1 

PLOT 205 

KNsO 

PLOT 210 

IN»0 

PLOT 215 

00 10 lal.NKT 

plot 220 

KN*KN*1 

PLOT 222 

IF d.EQ.DOOTO 20 

plot 225 

1F(KN.EQ.IDEl)QOTO 20 

Plot 2Jo 

READ(NPL)0UM 

PLOT 235 

GOTO 10 

PLOT 240 

INalNtl 

Plot 245 

READ (NPL) (PLO(IN.J) .U=1.H) 

PLOT 250 

5?t?N,EQ.NEND)GOTO 30 

Pb8i m 

continue 

plot 265 


20 


30 CALL PLOT! (PLO(l*D »PL0(1»2) f INtIfO,»,2»XNAMEl»YNAMXtPT»i»l.l,200)PLOT 280 

CALL PLOTl (PLO(l»n ♦PL0<1»3) . INt 1 » 0 , , , 2» XNAMEl # Y naM 2.PT t 1 t 1 ♦ I , 200 ) PLOT 290 
call PLOT! (PLO(l.l) .PL0(1»4) , I N. 1 » 0. » ,2 ♦ XNAME 1 . YNAM3 »PT ♦ 1 ♦ I ♦ 1 ,200 ) PLOT 300 
CALL PLOTi (PLOd.l) ,PL0(1.5) , IN. 1 . 0 . » ,2. XNAMEl ♦ YNAM4 ,PT , 1 , 1 , 1 , 200 ) PLOT 310 
CALL plot l(PLO( 1.1) .PL0(1,6) . I N. 1 . 0 , . ,2. XNAMEl . YN aMS ,PT . 1 . 1 . 1 , 200 ) PLOT 320 
CALL PLOtl (PLO(l.l) .PL0(1.9) , IN. 1 , 0 ... 2. XNAMEl . YNAM7 .PT , 1 . 1 , 1 , 200 ) PLOT 3 JO 
call PLOTl (PLO(l.l) .PLO(l.lO) .IN.lt 0,.,2. XNAMEl. YNAM8.PT.1. 1.1. 200PLOT 340 
$) PLOT 3b0 

CALL PLOTI (PLOd.l ) .PLOntll).IN,l»0,,.2,XNAMEl.YNAM9,PT.l.l,l,200PLOT 3&0 
S) plot 3/0 

CALL PLOTI (PLOd.T) .PLOdtB) . IN. 1 .XS, XO. XNAME2 . YNAM6.PT . 1 . 1 . 1 , 200 ) PLOT 380 
return PLOT 390 

END Plot 40o 
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SAMPLE PLOTS 





000-01 


8.000-01 


1 . 800+00 


1 . 600+00 


8.00 


VS SEC 

80DE74 LARGE AKPLITUDE SLOSH SltOjLATION 



)1 8 . 000-01 1 , 200+00 1 . 600+00 2.00 

S SEC 

LARGE A^■.PLITUDE SLOSH SIMULATION 




II 8.000-01 i.aoo+00 1.600+00 a. 00 

S SEC 

74 LARGE AMPLITUDE SLOSH SIMULATION 





II 8.000-01 1 .800+00 1.600+00 8.00 

S SEC 

74 LARGE AMPLITUDE SLOSH SIMULATION 
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)1 8.000-01 1.800+00 1.600+00 8.00 
S SEC 

74 LARGE AMPLITUDE SLOSH SIMULATION 
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C-1 


APPENDIX C - TEST AND ANALYTICAL RESULTS 


This appendix presents (Figures C-1 through C-22) the re- 
duced test data for all 22 test cases as well as the analytical 
prediction for each case assuming no viscous dissipative force 
(’?= 0, M= 0). In general, an update criteria (CRIT) of T’L was 
used, however, some cases were run with a value of 0.57o. It 
has been found that smaller values of CRIT can drastically im- 
prove the quality of the analytical prediction for some test 
configurations. This improvement can be seen by comparing 
Figures C- 17 and C-23. These figures both present data com- 
parisons for test 17 (50% fill, 45° tank inclination). The 
analytic results shown in Figure C-23 were run with CRIT = 27o, 
while those in Figure C- 17 were run with CRIT = 0.57o. The more 
lax update criteria (Figure C-23) allows excessive deviation 
from the constraint surface yielding more pronounced discon- 
tinuities at updates. 

Table C-1 summarizes the test conditions for each case. 
Values of axial acceleration (AZI) were calculated based on 
test time, drop capsule and drag shield masses, and drop cap- 
sule travel distance. Lateral accelerations (AYI) were scaled 
from high speed photographs taken during the testing. 

Table C-2 presents a qualitative evaluation of each test 
case including an assessment of the worth of the measured 
forces. Test eccentricities are also noted. 
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Figure C-17. Test 17; 50% Fill; ex • 45’; A = ,09g; CRIT ■ 0.5% 
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Test Data 

I . 500*00 
1 .000+00 
5.000-01 
0. 

-5.000-01 
-1.000*00 
-I .500*00 
-2.000*00 
-2.500*00 
-3.000*00 
-3.500*00 


0. 4.000-01 Q.OOO-Ql 1.200*00 l.G00*00 2.00 

t(sec) 



0. 4.000-01 □. 000-01 1.200*00 l.G00*00 2.00 

t{sec) 

1 .400*00 
I .200*00 
I .000*00 
0.000-01 
G. 000-01 
r 4.000-01 
2.000-01 
-3.553-15 
-2.000-01 
-4.000-01 
-6.000-01 

0. 4.000-01 0.000-01 1.200*00 1.600*00 

tisec) 

Figure C-23. Test 17; 50% Fill; ex = 45°; A = .09g; CRIT ■ 2.0% 
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TABLE C-1. AXIAL AND LATERAL ACCELERATION TABLE 



ex 


Test # 

(deg) 

7o Fill 



(in/sec^) (in/sec*^) 


0.0 

0.055 

0.50 

0.51 

1.10 

1.11 

2.045 
0.0 
0.06 
0.80 
0.81 

1.05 

1.06 
2.03 
0.0 
0.06 
0.60 
0.61 

1.45 

1.46 
2.03 
0.0 


386.04 

- 16.49 

- 16.49 

- 16.49 

- 16.49 

- 16.49 

- 16.49 

386.04 

- 16.79 

- 16.79 

- 16.79 

- 16.79 

- 16.79 

- 16.79 

386.04 

- 16.79 

- 16.79 

- 16.79 

- 16.79 

- 16.79 

- 16.79 

386.04 

- 15.83 

- 15.83 


4.0 

4.0 

4.0 
11.66 
11.66 

2.00 
2.00 
4.75 
4.75 
4.75 

- 1.66 

- 1.66 

3.80 

3.80 

2.81 
2.81 
2.81 
1.00 
1.00 
3.12 
3.12 
6.0 
6.0 


0.07 

2.08 


6.0 







C-26 



0.065 


2.075 


386.04 

- 16.40 

- 16.40 

- 16.40 

- 16.40 

- 16.40 

- 16.40 

- 16.40 

- 16.40 

386.04 

- 15.92 

- 15.92 

- 15.92 

- 15.92 

- 15.92 

- 15.92 

386.04 

- 16.69 


2.035 


0.065 


- 16.69 

386.04 

- 16.59 


- 16.59 

386.04 

- 15.92 

- 15.92 

- 15.92 

- 15.92 

- 15.92 


2.075 


- 15.92 
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TABLE C-1 (cont.) 



ex 


■asiM 

AZI 

AYI 

Test # 

(deg) 

7o Fill 

mSM 

(in/ sec^) 

(in/ sec^) 



2.9 

2.9 

2.9 

0.0 

0.0 

3.0 

3.0 

3.0 
5.71 
5.71 

1.87 

1.87 

1.87 

1.87 

1.87 
0.0 
0.0 

7.0 

7.0 
0.0 
0.0 

4.0 

4.0 

8.33 

8.33 

8.33 
0.62 
0.62 
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TABLE C-1 (cont.) 


Test # 

0X 

(deg) 

7o Fill 

Time 

(sec) 

AZI 

(in/ sec^) 

AYI 

(in/sec^) 

14 

0.0 

50. 

0.0 

386.04 

6.66 




0.065 

-34.18 

6 . 66 




1,62 

-34.18 

6 . 66 

15 

0.0 

75. 

0.0 

386.04 

1.90 




0.06 

-26.12 

1.90 




0.25 

-26.12 

1.90 




0.26 

-26.12 

3.50 




0.70 

-26.12 

3.50 




0.71 

-26.12 

-4.00 




0.93 

-26.12 

-4.00 




0.94 

-26.12 

-3.00 




1.22 

-26.12 

-3.00 




1.23 

-26.12 

0.0 




1.81 

-26.12 

0.0 

16 

45. 

25. 

0.0 

386.04 

8.33 




0.10 

-35.20 

8.33 




1.60 

-35.20 

8.33 

17 

45. 

50. 

0.0 

386,04 

6.25 




0.065 

-28.18 

6.25 




0.29 

-28.18 

6.25 




0.30 

-28.18 

0.0 




1.30 

-28.18 

0.0 




1.31 

-28.18 

7.0 




1.755 

-28.18 

7.0 

18 

45. 

75. 

0.0 

386.04 

6.0 




0.065 

-24.903 

6.0 




0.65 

-24.903 

6.0 




0.66 

-24.903 

-1.0 




0.92 

-24.903 

-1.0 
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TABLE C-1 (cont.) 


7, Fill 


Time 

(sec) 


AZI 

(in/ sec^) 


AYI 

(in/ sec^) 





TABLE C-2. QUALITATIVE EVALUATION OF DROP TEST RESULTS 


C-30 



TABLE C-2 (cont. 


C-31 








